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ABSTRACT 


IMPLEMENTATION OF BALDWIN-BARTH TURBULENCE MODEL 
INTO A TIME-ACCURATE CODE FOR UNSTEADY FLOWS 

Scott L. Low 
June 1993 

The Baldwin-Barth turbulence model was implemented into Zeta, a time-accurate, 
zonal, integro-differential code for incompressible laminar and turbulent flows. The 
implementation procedure patterned after the model subroutine in ARC2D. The results of 
ZETA with the Baldwin-Barth turbulence model were compared with experimental data, 
with ZETA using Baldwin-Lomax model, and with ARC2D using the Baldwin-Barth 
model. The Baldwin-Barth model subroutine was tested by inputting an ARC2D velocity 
solution of an NACA-0012 airfoil at Re = 3.91X10 6 and a = 5°. The resultant turbulent 
viscosity and Reynolds stresses compared favorably with the original data. For the same 
grid having grid points inside the laminar sublayer, which is necessary due to the one- 
equation nature of the model, ZETA however predicts early separation. It was found that 
the current ZETA has problem with such a fine grid. Further work is in progress to solve 
this problem. 
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NOMENCLATURE 


A Cross-sectional area of an airfoil. 

A + Constant equal to 26 in the turbulence models. 
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Bn Components making up By where n is an index. 

B s Internal boundary. 

Boo External boundary. 

bn Velocity Fourier coefficients where n is an index. 

Qj Upper elements of the tridiagonal matrix of the vorticity transport equation. 
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c n Velocity Fourier coefficients where n is an index, 
ci, C 2 Grid stretching parameters. 

Dij RHS matrix elements of the discretized vorticity transport equation. 
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Di J>2 Damping functions in the Baldwin-Barth model, 
dn Velocity Fourier coefficients where n is an index. 

E Extent of the grid. 

F(y) A function in the Baldwin-Lomax model. 

Fkieb Klebanoff intermittence function. 

jk A Damping function in the Baldwin-Barth model where k is an index or variable. 
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Transformation factor. 


JR Radial index of the grid. 

L, AL Airfoil chordlength. 

N^r) A M-type matrix operator in the discretized Baldwin-Barth model equation. 

P Production of k in the Baldwin-Barth turbulence model. 

R Flow region or radius of a unit circle. 
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Uoo Freestream velocity. 

u + Scaled velocity. 
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u’.v' Cartesian velocity components in the computational plane and inertial system. 
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v Velocity vector in the physical plane and inertial coordinate system, 

v Velocity vector in the physical plane and rotating coordinate system. 
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V P ,V 4> Cylindrical velocity components in the computational plane and inertial system. 
V P ,V 0 Cylindrical velocity components in the computational plane and rotating system. 
x,y Cartesian coordinates of the grid. 
y + Length scale. 


xiv 


z 


Physical plane: Z = x + iy, Z= re‘9 . 
a Airfoil angle of attack. 

a n Vorticity Fourier coefficients where n is an index. 

tta, a a The first coefficient s of the advective terms in the Baldwin-Barth model. 

a d, a d The first coefficient s of the diffusive terms in the Baldwin-Barth model. 

Vorticity Fourier coefficients where n is an index, 
ft*, P* The second coefficient s of the advective terms in the Baldwin-Barth model, 
ftd.ftd The second coefficient s of the diffusive terms in the Baldwin-Barth model. 
Am Increment of m where m is a dummy variable. 

5 Boundary layer thickness, 

e Dissipation terms in the turbulence models. 

”1* The third coefficient s of the advective terms in the Baldwin-Barth model. 
Td,? d The third coefficient s of the diffusive terms in the Baldwin-Barth model. 
k Von Karman constant or Production terms in the turbulence models, 
v Kinematic viscosity. 

v t> v e Turbulent and effective viscosities. 
v ti Inner turbulent viscosity in the Baldwin-Lomax model. 

v to Outer turbulent viscosity in the Baldwin-Lomax model. 

P.<t> Computational radial and angular coordinates. 

ct Physical plane origin shift in the transformation. 

? Computational plane: S = £ + C = P^. 

Q Body or airfoil angular speed, positive counter-clockwise. 

® Scalar and vector vorticities. 

General computational coordinates. 

t 

V Stream function in either plane and the inertial coordinate system. 

V Stream function in either plane and the rotating coordinate system. 


XV 


CHAPTER 1 
INTRODUCTION 


1 . 1 Motivation 

For years .unsteady flow has been the subject of numerous and continuing studies, 
in the field of helicopter or rotorcraft aerodynamics. The performance of a helicopter such 
as the UH-60A Black Hawks (Fig. 1.1) depends greatly on the aerodynamic lift of the 
rotor. Since the angle of attack of any rotor blade section oscillates as it travels through the 
rotor plane (Fig. 1.2), the resulting flow over the blade has complex and periodically 
changing characteristics. Understanding and accurately predicting the unsteady 
aerodynamics of flows over airfoils are, therefore, critical for designing new rotor blade 
airfoils and improving helicopter performance. 

1 .2 Previous Work 

Much experimental work has been done and many numerical flow solvers have 
been developed to study unsteady flows. Experimental data from unsteady flow 
investigations on geometries such as rotor system [1] and advanced airfoil sections [2-4] 
are readily available for review and comparison. Reference [2] contains the two 
dimensional dynamic stall characteristics of eight arifoils in sinusoidal pitch oscillations 
over a wide range of unsteady flow conditions. Numerical codes ranging from panel 
methods, coupling between time-dependent inviscid panel method and an unsteady 
boundary layer code [5], a full-potential code [6], an Euler code [7], a zonal integro- 
differential method [8] and other finite difference codes based on the Navier-Stokes 
equations, have produced results of reasonable agreement with experimental results. The 
code in Ref. [6] predicted comparable surface pressures for helicopter rotors and the Euler 
code of Reference [7] was used to investigate the rotor blade-vortex interactions. 
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2 

One recent study using the zonal integro-differential method for a two-element 
airfoil [9] and other single airfoils has demonstrated that this procedure is effective in 
treating general viscous flows, even with large separation regions. The code, ZETA [10] 
utilizes an integral representation of the velocity vector [11], a velocity-vorticity formulation 
of the Navier-Stokes Equations, and a Fourier series expansion. Using an integral 
representation of velocity, the flow computation may be confined to viscous zones. For 
turbulent flows, the Baldwin-Lomax [12] algebraic model is used to determine the eddy 
viscosity. A typical grid has 80x50 grid points, which is coarse compared to prevailing grid 
requirements for Navier-Stokes solutions. The general performance of this numerical 
procedure is satisfactory. However, further evaluations and refinements are required prior 
to using it for designing high-lift rotor blade airfoils. 

1.3 Present Works 

The Baldwin-Barth turbulence model [13] is a self-consistent one-equation model 
that does not require an algebraic length scale. It is derived from a simplified form of the 
standard k-e model equations [13]. This robust model was found to give significantly 
better results than the algebraic Baldwin-Lomax model in a recent study of four popular 
turbulence models [14]. This and other good results obtained with this model prompted the 
present investigation. The objectives are to implement the Baldwin-Barth turbulence model, 
and examine the effect of finer grids and Reynolds numberfRe) on ZETA. No known time- 
accurate case has be run with the Baldwin-Barth model. The computational results will be 
compared with that of ARC2D and with the experimental results of McAlister, et al [2]. 

This report presents the numerical procedures used and the results. Chapter 2 and 3 
describe the governing equations, grids, and boundary conditions. The code's numerical 
formulation is explained in Chapter 4, and the two turbulence models are presented in 
Chapter 5. Numerical results for the NACA-0012 airfoil and comparisons are presented in 
Chapter 6. Chapter 7 contains major conclusions and recommendations for future study. 



CHAPTER 2 


GOVERNING EQUATIONS 


2 . 1 Equations of Motion 

The behavior of a viscous, incompressible, turbulent flow is described by the 
continuity equation and the vorticity transport equation in the code. The vorticity transport 
equation is derived by taking the curl of each term of the incompressible Navier-Stokes 
equations in its familiar pressure-velocity form. Since the flow is incompressible, no 
energy equation is required. With body forces such as gravity and heat transfer are 
negligible, these equations can be expressed as 

Vv*= 0 (2.1) 


^=-(v'V)a+(w •v)? + V 2 (veG))+ S w 
dt 


( 2 . 2 ) 


where v denotes the velocity and the vorticity vectors: 

cd = Vxv (2.3) 

The effective viscosity v e is composed of both eddy viscosity v t and kinematic viscosity v: 

v e = v t + v 

The physical processes of convection, stretching and rotation, and diffusion of vorticity are 
represented by the right hand side terms of Equation 2.2 respectively. For a two- 
dimensional flow, the scalar source term S© reduces to 


I a 2 i 

du\ c 

2 ( dvl 

a 2 1 

av\ 

a 2 / 


2 -H 
dx 2 * 

i v %ra; 

^ Ve ax) 

+ dxdyl 

Ve ay) 

dxdy\ 

x, 

e a x ) 


The Sco term is negligible relative to the other terms of the equation. 

These equations are appropriate for both external and internal flow, but external 
flow will be emphasized. Equation 2.1 to Equation 2.3 are also three-dimensional but their 
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two-dimensional forms will be applied in ZETA. The corresponding velocity field is 
computed from the vorticity field with Equations 2. 1 and 2.2. Along with initial and 
boundary velocity conditions, these equations will uniquely determine the time-dependent 
flow field of an incompressible fluid. 

2.1.1 Kinematics and Kinetics of the Flow 

The equations of motions are divided into a kinematic aspect and a kinetic aspect. 
The kinematic aspect of the problem relates the velocity field to the vorticity field at any 
instant of time. For a given vorticity field, the velocity field can be uniquely determined 
[15]. This aspect of the flow is governed by Equations 2.1 and 2.3 which are linear and 
elliptic. The solution of these equations will require prescribed boundary conditions about 
the flow field. Since the flow field is known only at infinity, the entire flow field must be 
included in the solution procedure. 

With the vorticity- velocity formulation, Equations 2.1 and 2.3 may be reformulated 
as an integral representation for velocity vector at time t [1 1] 

v(r,t) = - 1 C 0 o xV o PdR o + <j> [(v 0 • n 0 )-(v 0 xn 0 )x]v 0 PdB 0 

A Jb (2.4) 

where B includes the internal boundary B s and external boundary B„ of the region R. n is 

the unit normal on B facing away from the region R. The subscript 'o' denotes that the 

operators and variables are in the viscous region, and P is the fundamental solution of 

Poisson's equation. 

PW ° )= 

With an integral representation, the integral over the fluid domain, R, does not need to 
include the inviscid region since C0o is zero there. The solution of velocity is then confined 
to the viscous region. 

The kinetic aspect of the problem deals with the change of the vorticity field with 


5 

time and it is described by Equation 2.2. This equation is nonlinear, parabolic in time, and 

elliptic in space. The equation describes the transport, not the generation or depletion, of 

— ♦ 

vorticity. Again, the solution of vorticity may be confined to the viscous region since C0o is 
zero in the inviscid region of the flow. This is the distinct feature of the numerical method. 
Taking advantage of this feature, the flow field is divided into three zones: an inviscid zone 
constituting the majority of the field, an attached viscous zone, and a detached viscous zone 
(Fig. 2.1). The detached viscous zone may include the wake, starting vortex assembly and 
separated regions of the airfoil. 

2.2 Coordinate Transformation 

The grid generation procedure employs a modified Joukowski transformation. The 
geometry and governing equations in Cartesian coordinates(x,y) are transformed to a 
generalized, body-conforming, curvilinear coordinate systemfP,^) [16]. The grid points 
have a one-to-one correspondence with the physical points. Unlike a conventional grid 
generator, this procedure works backward in that a specified computational grid is first 
constructed and then comformally mapped into the physical plane. The computational grid 
is composed of concentric circles and radial lines about a unit circle which represents the 
airfoil (Fig. 2.2). The grid transformation used is 

Ze ia = C + Y + c ~— + <* 

C +Y (2.5) 

where 

Z - physical plane 
= x + iy 
= re'6 

S = computational plane 
= £ + iri 
- pe^ 
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a = airfoil angle of attack 

c,Y = airfoil parameters 

ct = physical plane origin shift 
r,0 = physical radial and angular coordinates 

p,<{> = computational radial and angular coordinates 

The airfoil parameters are 

C = (£ + V 1 -T| 2 )( 1-5) (2.6) 

Y=S + ir l (2.7) 

The parameter c is the numerical chord length of the airfoil. Y is a complex translation of the 
origin to the center of the unit circle in the computational plane. 5 is a real number close to 
zero which prevents the transformation from becoming singular. It also specifies the 
curvature of the trailing edge. ^ will always be negative and "H will be zero for symmetric 
airfoils [15]. 

Based on the above transformation, the metrics and scale factor, H, of the 
transformations are defined as 


dx _ dy 
dx _ dy 

an~"a^ 


H = 


dZ 

ac 


an; 



( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 


The invariant of the transformation is the stream function. The integral representation for 
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the velocity vector is transformed into the computational plane by multiplying the vortxcity 

vector by the square of the scale factor H 2 , 

©(; = H 2 (0 (2.12) 

The vorticity transport equation is transformed by differentiating the physical coordinates 
with respect to the computational coordinates. The vorticity is computed in the rotating 
coordinate system so the grid parameters need not be recomputed as the airfoil angle of 
attack changes. The velocity is computed in a body-fixed inertial reference frame. When 
needed, the velocity values are transformed to the rotating reference frame. 

2.3 Transformed Equations 

The governing equations for the kinematic and kinetic aspects of the flow are 
presented in both the inertial and rotating coordinate systems. These equations are 
summarized in their differential forms. The two aspects of the flow problem are computed 
in different reference frames. Both coordinate systems are considered body-fixed. The 
transformed equations are presented in their respective frames of reference. In this section, 
a primed variable is associated with the inertial coordinate system while a non-primed 
variable is associated with the rotating coordinate system. TTie following list of variables 
are used: 

v Velocity vector in the physical plane and inertial coordinate system, 

v Velocity vector in the physical plane and rotating coordinate system. 

u',v' Cartesian velocity components in the physical plane and inertial system, 

u , v Cartesian velocity components in the physical plane and rotating system. 

u',v' Cartesian velocity components in the computational plane and inertial system, 

u ,v Cartesian velocity components in the computational plane and rotating system. 

v P> v <t> Cylindrical velocity components in the computational plane and inertial 
system. 

v P ,v <t» Cylindrical velocity components in the computational plane and rotating 
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system. 

V Stream function in either plane and the inertial coordinate system. 

V Stream function in either plane and the rotating coordinate system. 

Q Body or airfoil angular speed, positive counter-clockwise. 

r Position vector with respect to axis of rotation. 

The velocities defined in terms of the the stream functions are 


. d\J/ 

■ av 

u = 3~~ 

By 

Bx 

By 

d\|/ 

u '¥ 

V 3x 

. By 

v ._ av 

u "an 

as, 

d\|/ 

11 — - 

a^ 

an 

V ”'aiT 


(2.13) 


The governing equations of the code are presented in the inertial reference frame below. 
Kinematics: 


V v =0 

(2.14) 

Vxv = co 

(2.15) 

vb = Q x rg 

(2.16) 


vb is the velocity at any station on the body or airfoil surface. 
Kinetics: 


= -(v • V )© + V ^Ve©) 
at 


(2.17) 


Transforming to the computational plane, the kinematic governing equations become 


V C^ = ° (2.18) 

I I 

V;xv^ = 0); (2.19) 


The two dimensional vorticity in the computational plane is 


Numerically, the kinematic aspect of the flow will be solved using the integral 
representation of the velocity vector and a Fourier series expansion. The governing 
equations of the code are now presented in the rotating coordinate system. 
Kinematics: 


Kinetics: 


< 

II 

o 

(2.20) 

V x v = co 

(2.21) 

v = v - £2 x r 

(2.22) 

co = co - 2£2 

(2.23) 

& 

II 

o 

(2.24) 

(v- v)c 0 + V^VeG)) - 2Q 

(2.25) 


The transformation of the governing equation of the kinetic aspect of the flow will be 
presented in more details. With the above equation and Equation 2.2, the vorticity can be 
referenced to the inertial system. The vorticity vector has been replaced by its scalar value. 


^- = -(v- Vjto’ + V^vcto') 
dt 


(2.26) 


The vorticity is in the physical plane and the velocities are still referenced to the rotating 
coordinate system. The divergence of velocity and Laplace operators are transformed as 
described in reference [10]: 


7 ' V = h^' V; ) 


(2.27) 


V 2 =-i-V 
H 2 


2 


(2.28) 


The resulting vorticity transport equation in the computational plane and in its conservative 
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form is 


H2 ^- + V ; (v ; o)') = v ^ v e co) 

To obtain the numerical form of this equation, it is first non-dimensionalized. 


(2.29) 


x * = x 
R 

u*=-U- 

U- 

<o*=ak 


y R 
v* = JL 
U- 
* = tU- 


U- ' R 

R = 1 is the radius of the unit circle in the computational plane which represents the airfoil 
surface. U- = 1 is the velocity at infinity. The non-dimensionalized vorticity transport 
equation is 


(2.30) 


9f 

With L defined as the airfoil chord length, the Reynolds number Re is 

R U ~ L 

Ke_ v (2.31) 

The numerical kinematic viscosity v is therefore defined as L/Re since U- = 1 . Dropping 
the superscripts and writing in cylindrical coordinates, the vorticity transport equation to be 
discretized becomes 

■>9© 


1 9 [ P 9COy ] | 1 9 2 C0y 

P9* 2 


(p-c 2 )9s\(p-C2) 9s 


(2.32) 


where 


CDv 


= fe +v ‘)“ 


Re 

9 _ 9s 9 
9p 9p9s 


P = e (s+Cl)+C2 


The variable s, ci, and C2 will discussed in details in the next chapter. 

Finally, the kinetic part of the problem requires velocity values in the rotating 
reference frame. Since the velocity is solved in the inertial reference frame in the kinematic 


part of the problem, the velocity correlations in the computational plane are given here. 
Having obtained the cylindrical velocity components, v p and v <t>, in the inertial coordinate 
system, the Cartesian velocity components, u’ and v', are computed as 

u' = vj,cos (J> - v^sin <|> 


v' = VpSin <}> + v^cos <|> 


(2.33) 

(2.34) 


The Cartesian velocity components are then transformed to the rotating coordinate system. 


, 3x„ 3y~ 

u = u + 5iT^ + 5i^ 


(2.35) 


v = v' - — Q x - — £2y 

^ ^ (2.36) 

Transforming back to the cylindrical coordinates, the velocity components in the rotating 
coordinate system become 


v p = ucos <(> + vsin <|> 
v^ = vcos <J> - usin <f> 


(2.37) 

(2.38) 


2.4 Aerodynamic Loads 

Use of the vorticity-velocity form of the Navier-Stokes equation has the 
disadvantage of not having the pressure distribution computed from the equations of 
motion. Thus, special attention is required . The total head, h, of an incompressible flow 
can be derived from the pressure- velocity form of the Navier-Stokes equations [10]. An 
integral representation of the total head, h, for a two-dimensional external flow is 



(v o x0)o)(r o -f) 

M 2 


dRo 



h 0 (r 0 -r) n 0 
^o-ll 2 


dB c 


co 0 (f 0 -r)s 0 
^o-ll 2 


dB 0 + h,. 


1 2 


(2.39) 


where B s is the internal boundary, and n and s are the unit normal and tangent vectors, 
respectively. The pressure coefficient, Cp, is calculated as 


Cp(p=l,<t» = h(p=l,<» - h» + 1 (2.40) 

where P and <1> are the radial and angular coordinates of the computational plane. Having 
determined the pressure coefficient, the components of aerodynamic force coefficients on 
the body(P = 1), in the computational plane, are computed as follow: 

f 2n 

C Np = -f Cp«)))^d<t> 

^ d<J> 


r2n 


c N f =-^-| CD(l,<j>)^d<|> 

Re I 0<J) 


(2.42) 


rln 

C Tp = f| Cp(l,0)^d<D 

L * a<i> 


(2.43) 


fin 


Or = -2- 
Tf Re 


C0(l,<|>)^d<|> 

a<t> 


(2.44) 


Cn = Cn p + Cn, 


(2.45) 


The subscripts N and T denote the normal and tangential components, and p and f the 
pressure and friction components. The lift and drag coefficients are 

Cl = Cncos a - Opsin a (2.47) 

Cd = Cnsui a + Ctcos a (2.48) 

where a is the airfoil angle of attack. The moments are computed with respect to the 
quarter-chord and are positive counterclockwise since ZETA uses a left-handed coordinate 
system. The moment coefficients are 


Cm p = 
Cm, = 


L 2 


Rc-L 


Cp(<t>) 

' 9x dy" 
x — + y— 

d<(> 


3<j) 


( ®( 1 ><t>)f x— - y— 1 

A 

L 

9<t>j 


d(J> 


(2.49) 


(2.50) 


Cm = Cm p + Cm, 


(2.51) 


All of the above coefficients are normalized with respect to the airfoil chord length L. 


CHAPTER 3 


GRID AND BOUNDARY CONDITIONS 
3 . 1 Grid Generation 

The grid generation procedure as mentioned in Section 2.2 uses a modified 
Joukowski transformation. Though a numerical transformation is also available, only the 
well tested Joukowski transformation will be discussed. This transformation conformally 
maps the computational plane about the specified airfoil in the physical plane. The grid is 
composed of concentric circles and radial lines in cylindrical coordinates about a unit circle 
which is the transformed airfoil surface. The grid generator constructs the physical and 
computational grids, and computes the parameters of the grid transformation and all the 
variables which depend only on the specified grid. 

The radial lines in the computational plane are evenly spaced about the unit circle. 
When conformally mapped into the physical plane, the radial lines arc more concentrated 
near the leading and trailing edges (Fig. 3.1). The concentric circles in the computational 
plane are, however, stretched in the radial direction. With conformal mapping, the 
stretching will give better resolution near the leading and trailing edges. The grid generator, 
presently, does not have any adaptive feature like clustering, but the conformal 
transformation accounts for it by concentrating lines near the surface and the leading and 
trailing edges. 

The stretching formula in the computational plane is 

p = e s + Cl + C2 (3.1) 

where s = (n - l)As 

n is the index of the radial grid line and s is the distance from the unit circle. The radial grid 
extent E in the computational plane and the desired spacing between the first two azimuthal 
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lines Asi must be specified in order to determine the unknown stretching parameters As, ci 
and C 2 - A minimal non-dimensional grid extent of 18 that corresponds to approximately 5 
chordlengths is recommended. Typically, the boundary layer thickness and the desired 
number of grid lines within the boundary layer will determine the first grid spacing Asx. 

For some applications, the Asi will have to be inside the laminar sublayer. The laminar 
boundary layer thickness is [17] 

g = 5AL 

V1T (3.2) 

where AL is the non-dimensional airfoil chordlength. The turbulent boundary layer 
thickness is 

g .37AL 

Re V5 (3.3) 

Knowing the grid extent E, the radial grid dimension JR, and the first grid spacing Asi, 
along with defining the unit circle will give a system of three equations to determine the 
three unknown parameters. An area ratio of less than 1.2 was maintained and experience 
has shown that As should be less than 0.15 for an effective grid [10]. The system of 
equations is 

e Cl + C 2 = 1 .0 

c As-k:i + C2 = 1.0 + Asi 

gJRAs+ci + c 2 = E (3.4) 


3.2 Boundary Conditions 

For a given airfoil geometry, the combination of initial and boundary conditions 
distinguishes the flow patterns. The airfoil and the surrounding fluid are initially at rest. 
Immediately after time t = 0, the airfoil instantaneously translates at a velocity of -v~. The 
flow is potential at this instant. A sheet of vorticity acts as a velocity discontinuity between 


the airfoil surface and the undisturbed surrounding fluid. This vorticity sheet will diffuse 
and convect away with time. 

For the solution of the vorticity and velocity vectors, four boundary conditions are 
required. These are the vorticity and velocity boundary conditions on the flow field s 
internal and external boundaries. Each of these is a Dirichlet type boundary condition. The 
internal boundary is the airfoil surface and the external boundary is taken to be at infinity. 

The velocity on the external boundary is zero since the fluid far from the airfoil is at 
rest. The external boundary condition is satisfied exactly. The v- term in Equation 2.4 
accounts for the relative velocity between the airfoil and the fluid. The velocities on the 
internal boundary are known through the no slip condition and the prescribed airfoil 
forward speed, angle of attack, and oscillatory motion. The integral over the internal 
boundary in Equation 2.4 will generally be non-zero since the airfoil is free to rotate about 
some body-fixed origin. This boundary integral requires both normal and tangential 
velocity components, but one is sufficient to uniquely determine the incompressible 
velocity field in R. 

At the external boundary, the vorticity is simply zero. The vorticity external 
boundary can be located anywhere inside the inviscid zone. Computation of vorticity is 
confined to just inside the inviscid zone. A zero vorticity gradient normal to the boundary is 
applied when cutting through the vortical wake is required. The internal boundary 
vorticities are accurately computed with Equation 2.4 which is part of the kinematic aspect 
of the problem. Knowing either the normal or tangential velocity component on the airfoil 
surface, the near-surface vorticity co can be computed uniquely [18]. However, the same 
reference shows that using the normal velocity component may introduce numerical 
difficulties, so the tangential component is used. 


CHAPTER 4 


NUMERICAL METHOD 

The numerical formulation of the governing equations is an integro-differential and 
zonal methods. Applying the zonal procedure does not require modification of the 
governing equations. As previously stated, the kinetic and the kinematic aspect of the flow 
are treated separately through the integro-differential approach.The great advantage of using 
this approach is that it permits the flow solution to be confined to the viscous regions. The 
kinematic aspect or the velocity vector equation utilizes the integral approach while the 
kinetic aspect or the vorticity transport equation utilizes the differential approach. There are 
three major steps in the computation loop. The loop starts with one kinetic part and 
followed by two kinematic parts. It may be summarized as [7] 

1) Solving the vorticity transport equation, the interior vorticity values at 
the new time level are computed with the known vorticity and velocity 
values at the previous time level. 

2) Using the newly computed interior vorticity values, new boundary 
vorticity values are established. 

3) A new velocity field is computed with the new vorticity field. 

The velocity is computed at grid points and the vorticity is computed at half points in the 
radial direction. Details of the numerical methods used for the two aspects of the problem 
will be discussed. 

4.1 Kinematic Aspect Of The Flow 

The integral representation for the velocity vector is 
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211 P.4P 

y r 



dB 0 + Vo, 


(4.1) 

The boundary conditions and two-dimensional restrictions have been imposed on this 
equation. Since the computational grid is in polar coordinates, the velocity components may 
be expressed as 


vr(r,0) = -L 


0)oroSin(e o -8) 


dRo 


I r§+r 2 -2r o rcos(0 o -0) 

' R 
*2% 

v r Ir 0 cos(8o-8)-r] rf0 
^ 71 1 r^+r 2 -2roTcos(0 o -0) 


ve/oSin (e o -8) 

I r£fr 2 -2r o rcos(0 o -0) 


Tod0 o + V 0 


(4.2) 


v,( r> 9)-^| 4r^(e. : 9)-r] 

I r£+r 2 -2r 0 rcos (0 O -0 ) 


r R 

,2 k 




VroToSin (Qq-Q) 

2 K J r?+r 2 -2rorcos(0 o -0) 

f | v4r 0 cos(9 0 -9).r] ^ 

I rg+r 2 -2r o rcos(0 o -0) 


(4.3) 


where 
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Voo, = v„cos(a)cos(0) + v oo sin(a)sin(0) 

= v o „sin(a)cos(0) + v oo cos(a)sin(0) 

and to is again replaced by it sole component co. As in Equation 4.1, the last two integrals 
of Equation 4.2 and 4.3 are over the interior boundary B s which is the airfoil surface. 

These two equations are also valid in the computational plane if is replace by which 
is transformed as 

= (4.4) 

This transformation was discussed in the last chapter. In the computational plane, the 
position components r and 0 are replaced by P and <t>. The angle of attack a is positive 
nose-up as before. 


4.1.1 Fourier Series Expansion 

Using conventional methods to compute the integrals for the velocity components 
will be inefficient numerically. If the finite Fourier series expansion is used for the 
integrals, the vorticity and the velocity components can be evaluated explicidy. The 
modified equations are 


®^P><1>) = ^4^ + X (otntpJcosCn^+PnfpJsin^)) 

L n=l 

+ ^Plcos(N<))) 

4r 

v p (p,4>) = Np- + X (a n (p)cos(n<(>)+bn(p)sin(n<|))) 

Z n=l 

+ ^^COS(N<{>) 

v$(p,<\>) = ^4^- + X ( c n(P ) c °s(n<l))+d n (p )sin(n<l>)) 

L n=l 

+ ^^cos(N<j>) 


(4.5) 


(4.6) 


(4.7) 
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Each of the Fourier coefficients is a function of P only, assuming the coefficients are 
constant along any azimuthal grid line. There are 2N = 1-1 terms in each series. I is the 
total number of discrete points in the tangential direction about the grid. Since an "O" grid 
is used, the first and last points overlap. Knowing ( 0 5 at all grid points, in particular, at 2N 
points along an azimuthal line, the Fourier coefficients can be determined as 

2N-1 

a k = 7 f X a)£(<|>p)cos(k<l>p) k = 0,1,...,N 

p=o (4.8) 


2N-1 

Pk = iX co ; «|>p)sin(kV k = 0,1,...,N-1 

p=° (4.9) 

Substituting Equations 4.5 through 4.7 into Equations 4.2 and 4.3 in the computational 

plane and evaluating the integrals, the following relations between the known vorticity 

Fourier coefficients and the velocity Fourier coefficients are: [12] 


a<>(p) = 


ao(l) 

P 


(4.10) 


ai(p) = ^| Pi(po)(y) 2 dpo + ^J Pi(Po)dp c 
+ ^W[ai(l) + di(l)] + v 00 cos(a) 


(4.11) 


&n(P> = ^ PnCpo^^J dPo + 2 f Pn(po)|~|^ dp 0 
^ 0^1 2<n<N 

bi(p) = ai(po)jy) dp 0 - if ai(po)dp 0 

'^pj 2[ci(1) ' bl(1)] + v “ sin(a) 


(4.12) 


(4.13) 


b n (p) = 2 f °^ (po) {^rf <ip ° " 2^f ° tn(po) (^') n dp ° 

2 1 

J 1 W P 

2 £ n <N 

Apl 

(4.14) 

C 0 (p)=j a o(Po)|y)dpo + °p ) 

(4.15) 

ci(p) = ai(po)Jy)dpo -^1 ai(po)dpo 


+ ^-pcia) * biC 1 )] + v ~sin(a) 

(4.16) 

Cn(p)=^f dp ° * 2 f Otn(po) (^') n dpo 


W 1 

+ ilir , [ c " ( l)-b»a)] 2SnSN 

2\P/ 

(4.17) 

dl(p ) = 2 j ^ l(po) (^) P l(po)dpo 


+ x(— f[ai(l) + di(l)] - v^cosfa) 

API 

(4.18) 

d„<p) = 2"f -d Pn(Po)(yP<JPo 


J 1 »P 

+ ifirW 1 ) + d n(l)] 2 < n < N 

2\P/ 

(4.19) 


an(l), bn(l), Cn(l), and dn(l) are the known Fourier coefficients of the transformed 
velocities on the unit circle. By applying these Fourier coefficient equations on the 
tangential velocity component on the unit circle, the constraints on the vorticity Fourier 
coefficients are found. 


I ai(po)dp 0 = -ci(l)-bi(l) + 2v co sina 

)> (4.20) 
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Pi(Po)dp 0 = -di(l) + ai(l) - 2v eo cos a 


f 

l a n (p 0 )(J-) n l dp 0 = -c n (l)-b n (l) 

| P„(po)|^-) n ' 1 dpo = -d n (l) + a n (l) 


(4.21) 


(4.22) 


(4.23) 


The principle of conservation of total vorticity is employed to add the final constraint to the 
kinematic aspect of the flow. Using the finite Fourier series expansion of vorticity (Eq. 
4.5), the equation can be re-expressed for numerical application 



i 


C0 o dR o = -2A£2 


a 0 (po)Podpo = ~ Q 

7C 


(4.24) 


where A is the airfoil cross-sectional area and Q the airfoil angular speed. The velocity and 
the vorticity field are uniquely determined with these constraints. 


4.2 Kinetic Aspect of the Flow 

A combination of finite difference schemes is used to compute a numerical solution 
of the kinetic part of the problem in the computational plane. Discretization of the vorticity 
transport equation (Eq. 2.32) generates a set of finite difference equations(FDE) which are 
algebraic. The working variable is the vorticity as mentioned. As time advances, the 
employed time-accurate numerical method solves for the vorticity at the half grid points in 
the radial direction, based on known vorticity values at the previous time level. A first order 
accurate backward difference scheme is used on time, and a second order accurate central 
difference scheme in space is used on the diffusive terms. An upwind difference scheme is 
used on the convective terms. The differencing is second order accurate in the stream- wise 
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direction and first order accurate in the reverse direction. Collecting the terms of the FDE 
into a tridiagonal matrix for successive line under-relaxation along the radial line, the FDE 
becomes 

AijCiij- 1 ! + B.jWg 1 + GjWij-i = i,co? + i j) (4.25) 

The superscripts, k+1, indicate the new time level. The matrix elements for the laminar 
terms are 

Ay = A2 + A4 

= p,+m [v P R<o] — ~ p>+m -^ 

( Pj -c 2 )As (p r c 2 )As 2(pj+1/rC2) ( 4.26) 

Bij = Bl +B2 + B3 + B4 + B5 


where 
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+ 21VR, 
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Cij = C2 + C4 

Pj-l/2 


[VpL> 0 ]-- 


L/Re Pi-«g 


(Pj*c 2 )As L ' J (p r c 2 )As 2 ( p j' 1/rC2 ) 
Dij = Dl +D3 + D5 

= ^-pj(co^j) + (co^+i J^V^r<0] 

At A<|) 

- -^^-ijv^] + (trf + , j + crf-u) 


A<(> 


pjA<t> 

VpR = Vp u+1 
VpL = Vp y 


(4.27) 


(4.28) 


(4.29) 
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V <t>R ^ v <*!g +V !>u+i +V( t , i*u +V< l>t-*i>ij 
v <fL = 

The index j and J are for the vorticity and velocity grids respectively. The conditions in the 
brackets indicate when the term is included. The relationships between some of the 
variables are 

B2 = C2 
B3 = A2 
B4 = -A4 - C4 

B5 = 2 x D5 (4.30) 

When the flow at a grid point is considered turbulent, the diffusion terms, A4, B4, B5, C4, 
and D5, are multiplied by 



This modification is possible by utilizing the definition of ®v used in Equation 2.32. In the 
boundary layer, the discretized vorticity transport equation may be simplified to the 
boundary layer equation which is parabolic, and the solution is marched forward without 
iterating. With the zonal approach, the numerical procedure is made even more efficient. 

The numerical procedure used for the aerodynamic loads are very similar to that of 
the velocity vector equation in its integral form. Greater details are found in reference [10]. 


CHAPTER 5 


TURBULENT MODELS 

In examining the differential form of the incompressible turbulent, Reynolds- 
averaged Navier-Stokes momentum equation, there are apparent stress gradients due to 
transport of momentum by turbulent fluctuations and deformations attributed to 
fluctuations. To solve these equations by a finite-difference method, some closing 
assumptions have to be made about these apparent turbulent stresses. In 1877, Boussinesq 
suggested that the apparent turbulent shearing stresses might be related to the rate of mean 
strain through an apparent scalar turbulent or "eddy" viscosity [16]. Many turbulence 
models have been developed but all have limitations. The models range from algebraic to 
k— e formulations and the accuracy of these models' predictions varies. 

The present code employs the Baldwin-Lomax algebraic turbulence model [12]. 
This model requires small amount of computational time and has been used extensively. Its 
accuracy has been found to be comparable to more complex turbulence models [20]. The 
accuracy of the model will change depending on the flow conditions since it was calibrated 
and verified with experimental data. The present code using the Baldwin-Lomax turbulence 
model has been generally under-predicting the aerodynamic loads such as Cl. This differs 
from most turbulence models which over-predict the value of Cl- The model is used for 
attached and separated flow regions. Beside the inherent shortcomings of the zero-equation 
model [16], the calculation of the length scale will not be accurate for separated flow 
regions. To further develop the code and to improve upon its performance, the Baldwin- 
Barth turbulence model is to be added to the code. Since the major aim of this investigation 
is to implement the Baldwin-Barth one-equation turbulence model, this model will be 
presented in greater details. 
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5.1 The Bald win-Lomax Turbulence Model 

The Baldwin-Lomax turbulence model is a two-layer, zero-equation, or algebraic 
eddy viscosity model. This model is patterned after that of Cebeci [21] with modification 
that avoid the necessity for finding the edge of the boundary layer [12]. The eddy viscosity 


v t is computed as follow 



y — y crossover 
y crossover ^ y 


(5.1) 


where the subscripts i and o denote the inner and outer layers respectively. y C rossover is the 
smallest value of y at which the eddy viscosity values from both the inner and outer regions 
are equal. 

The inner eddy viscosity Vu is based on the Prandtl-Van Driest formulation. 

Vti = 1 M (5.2) 


where 


1 =ky[l- exp(-y + /A + )] (5.3) 

M is the vorticity magnitude and y + is the length scale. 

y v (5.4) 

where ux is the friction velocity and v is the kinematic viscosity. The friction velocity is 
used to scale the tangential velocity. The scaled tangential velocity is given as reference [22] 
which is the 'law of the wall'. 


in viscous layer 
u + = j4n y+ + C in inertial layer 


u + = y+ 
— 1 


(5.5) 


C is found by assuming that y+ = 10 is the matching point between the two layers. 

The eddy viscosity in the outer region, v to is given as 

Vto = KCcpFwakeFidetjfy) (5.6) 


with 
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Fwake “ ^l^(ymaxF max? C w kY max Udif/Fn 


(5.7) 


The value of Fmax is determined from this function: 

F(y) = yN [ 1 - exp(-y + /A + )] (5.8) 


and ymax is the corresponding y location. Fy e b is the Klebanoff intermittency function: 


Fkleb(y) = 



(5.9) 


2 

The quantity u dif is the difference between the maximum and the minimum total velocity in 
the profile. The constants used are 


A + = 26.0 C cp =1.6 Ckieb = 0.3 


C w k=10 k = 0.4 K = 0.0168 


Numerically, the model begins with the determination of friction velocity, u x , 
iteratively. Consequently, the inner and outer eddy viscosity, v t i and v t0 , is computed 
along the direction normal to the surface. The final eddy viscosity value is assigned 
according to the ycrossover location. 


5.2 The Baldwin-Barth Turbulence Model 

The originally proposed Baldwin-Barth turbulence model [13] is chosen over the 
recently modified formulation [23] because applications of the model showed that the 
original model formulation gave significantly better results [14]. The finite difference 
method is used to determine the turbulent Reynolds number Rj of the partial differential Rj 
equation. The turbulent Reynolds number is directly related to eddy viscosity v t . The 
derivation of the k-Rj and then the Rt model equation begins with a standard form of the 
k-e equations. 


_ v (v+^-Wk + P - e 
Dt \ oj 

(5.10) 

= V -Iv+^lVe + c - cJ^- 

Dt \ c E ) 1 k 2 k 

(5.11) 
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where the total or substantive derivative is 

J = | + V-V 

5t (5.12) 

and P is the production of K in the equation, 

P=v & t ®.ZvH 

\0Xj dXi/dXj 3 ldx k ; (5.13) 

By considering Rj and its differentials, 

p _ K 2 

T ve (5.14) 

dRr _ 2 dK . d£ 

Rt k e (5.15) 

a k-Rj equation is derived from the K-e equations with a valid simplification of the 
diffusion terms. 


= (2-c El )^Ip + (c Ei -2)k + |v+^-Jv 2 (vR t ) - £(Vv) V(vR T ) 


(5.16) 


where 


v t = Ch(vR t ) (5.17) 

and vRj is the appropriate field variable rather than Rj. To obtain the Rt equation, it will 
first require the rearrangement of Equation 5.14 

£ = J£ 2_ = {KijJC2£ 

vRt vR t (5.18) 

where K 2 = (K 1 +K 2) 2 and assigning a value to Kj without loss of generality, at large Rj- 


Kj = vRtP 

With the above two equations, a relation among k, ki, and k 2 is obtained. 

K = VvEiF(i+^) 

Substituting Equation 5.20 into Equation 5.16 and rearranging the terms leads to 


(5.19) 


(5.20) 
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3^ = (cgj - c a )W + (v + ^l|v 2 (vRt) - J-(vv,)- V(vR t ) 
- (2 - c J ** VvRtP - (2 - C t >2 

Ki + K2 


(5.21) 


Neglecting the last two terms in this equation, a self-consistent one-equation model is 
obtained [13]. This equation is valid over a major portion of the shear layer at sufficiently 
high Rj. For the model to be applicable in near-wall regions, the turbulent Reynolds 
number, Rt is re-expressed as 

R T = Rrf3(Rr) (5.22) 

where f3 is a damping function that approximates % ~ Rt at large Rt- Applying 
commonly used damping functions in the k-e models [13], the eddy viscosity is 

v t = vc^R T = vc^R T (5.23) 

and 


e= e -D .jg-.fcLLSf, 

VRt vf3Rj 

Kj is also applicable at small Rj when Rt is replaced with Rj in its definition. The 
resulting model equation for Rj is 

5^ = (c J 2 - c £l + (v + ^vIvRt) - ^vRt) 


(5.24) 


a E c e 


(5.25) 


Subjecting this equation to the thin shear layer approximation and further developments 
[13], the following list of functions are obtained to determine vRj. 


= ~^{ c e2 ’ c ei)^ 

O e 

V t = C^(vR T )DiD2 

Pt = pv t 

Di = 1 - exp(-y+/A + ) 
D 2 = 1 - exp(-y + /A 2 ) 


(5.26) 

(5.27) 

(5.28) 

(5.29) 


(5.30) 


Ur) =ct + (i-^)( K ^+ d,d 2 )(vi5^ 

+ t ^ |UI (^ x p(-yV A 2) D 2 + ^«p('Wai)oi)) 

The constants used for the preceding equations are 

k = 0.41 Ce. = 1.2 c E2 = 2.0 

c n = 0.09 A + = 26 A 2 = 10 


5.2.1 Numerical Formulation 

The numerical solution of the one-equation model is de-coupled from the flow 
equations. An implicit factored ADI solver for scalar equation is employed on a two- 
dimensional, logically rectangular mesh. This computational grid, which is different from 
the cylindrical computational grid of the main flow solver, is chosen to simplify and 
quicken the implementation process. The finite difference formulation of the model 
equation is patterned after that of ARC2D [24]. Having to use an "O" grid and computing 
time-accurately requires, respectively, periodicity and modifications for a time-accurate and 
integro-differential code for incompressible turbulent flows. 

Defining the solution vector R where Ry ~ vRifxij, y.j), the model equation with 
discretized advection and diffusion terms becomes a system of ordinary differential 
equations of the form [13] 


R, + m(r)r = DR 


(5.34) 


where M[R] is a M-type matrix operator representing the discretization of advection and 
diffusion terms and D is a diagonal matrix. A M-type matrix is diagonally dominant, has 
positive diagonal entries and negative off-diagonal entries, and has zero row sum. 


Since the source term of the one-equation model can be computed explicitly once 
the velocity gradient is known, only the advection and diffusion terms are discretized. 
Assuming v is constant, the model equation may be rewritten as 

V- V(vR t ) = 2 fv + JV(vR t ) - J-V-v,V(vR T ) + (c c J 2 - CeJVcJu/jSRr 
Ot \ a E l o E (5 35 ) 

where S is 

\dXj dxj/dxj (5 36 ) 


The advective terms is approximated by first-order accurate upwind differencing: 

V- v(vR t ) - aJRj+ij + paRij + yJRi-i j 

+ CC^Rij+i + PaRij "t" TfaRiJ-1 


where 


,x _ 1 ir 

7a = 

1 Uitj 
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Ax 
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Ay 


Ay 


(5.37) 


P![ = -(a a + 7 t) 

P? = -te+tf) 

The diffusive terms are approximated by second-order accurate central differencing and 
may be combined as 


where 


2 (v + ^(v^vRt) - V.^-V(vRt) - ctjR WJ + pX + iJRhj 

+ a^Rij+i + pdRij + TdRiJ-i 


(5.38) 


<*d = 

f d = 


Ax 2 

Ax 2 


W v+ -) - 

(-) 

n cj/jj 

'o Ji + l/2jL 

Wv + vt) - 

(-) 

n o'ij 




Pd = -(a5 + 75) Pd^^ + Tj) 

The matrix operator will not be a M-type matrix due to the coefficients of the diffusive 
terms if the grid resolution is poor. This condition is stricdy enforced in the algorithm. 


Typically, a grid wall or first grid spacing of y + = 3.5 is required since Rj is designed to 
behave linearly in the near- wall region for zero gradient boundary layers [13]. Therefore, 
grid of adequate resolution is critical to the model's performance. The recommended 
boundary conditions are 

1) Solid walls: Rj = 0. 

2) Inflow: Specify Rx to match experimental v t. 

3) Outflow: Extrapolate Rj from interior values. 

4) Freestream: Set to a small value < 1. 


A listing of the one-equation model is in Appendix A. 



CHAPTER 6 


RESULTS AND DISCUSSION 

Implementation and verification of the Baldwin-Barth one-equation turbulence 
model and evaluation of ZETA are described. ZETA is a zonal integro-differential code for 
incompressible, laminar and turbulent flows. The emphasis will be on the verification and 
evaluation processes since the field equations and numerical method for the implementation 
have been described in Section 5.2. The new constraints introduced by the Baldwin-Barth 
model and modifications of the ZETA require further evaluation of the code. The test cases 
were chosen for their simple flow characteristics and for correlation with the experimental 
results [2]. Numerical results are first compared with corresponding results of ZETA with 
the Baldwin-Lomax algebraic turbulence model and of ARC2D using the Baldwin-Barth 
model. The comparison will also utilize experimental results when applicable. 

Computations are performed on a NACA-0012 airfoil which has been extensively 
studied. Most cases have a Reynolds number Re of 3.91x106, an angle of attack a of 5.0 
degrees, and a fireestream Mach number M« of 0.301. This flow condition corresponds to 
an experimental case of Ref. [2], It was chosen to minimize compressibility effects so 
using the incompressible flow solver, ZETA, was justified. 

This flow case with simple attached flow features is selected so the newly 
implemented Baldwin-Barth subroutine can be easily tested. Although ZETA is capable of 
computing unsteady flow over oscillating airfoils, the computed cases are steady at a given 
angle of attack. The code is time-accurate and both laminar and turbulent computations 
were performed (Table I). All computations were done on the CRAY Y-MP 8/832 or C-90 
but the code can be run on VAX, IRIS, or SUN workstations. 

ZETA can only compute on the "O" grid which is created with a grid generator that 
is used exclusively for this code. The grid generator also provides the transformation 
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factors and other grid dependent variables for the flow solver. Creating a grid normally 
requires a CPU time of less than 1 .0 second. The computation begins when the airfoil 
starts impulsively from rest at a time immediately after t = 0. For most flows, the solutions 
are considered steady at a non-dimensional time of about 150. This time is equivalent to 
having the airfoil travel 42 chord-lengths through the fluid since the non-dimensional 
freestream velocity is 1 .0 and the chord length of an airfoil is 3.6190589 due to the 
conformal mapping scheme. The typical number of iterations to reach a non-dimensional 
time of 150 is about 7,600. This corresponds to a total CPU time of about 4.5 hours, or a 
real time of about 2 to 3 days. The convergence rates vary with grids, angles of attack, and 
Reynolds numbers. 

6. 1 Implementation of the Baldwin-Barth Turbulence Model 

The implementation of the Baldwin-Barth one-equation turbulence model is 
patterned after that of ARC2D [13]. Modification of the ZETA code was also necessary. 
The major change in ZETA involved the computation of physical velocity for the Baldwin- 
Barth subroutine. Presently, the model subroutine is capable of handling the periodicity of 
an "O" grid, the left-handed coordinate system, and the time-accurate aspects of ZETA. 
Since the model subroutine is de-coupled from the flow solver, a conventional rectangular 
computational grid was chosen for convenience. 

The numerical procedures in the model subroutine include computations of the 
metrics and Jacobians of transformation [16], a generalized distance function which gives 
the minimum distance to the solid wall [25], and the gradient of the mean velocity. The 
solution of the model equation utilizes an alternating-direction implicit(ADI) scheme in a 
half-stagger grid. The eddy viscosity is provided at half points in the radial direction. The 
"O" grid is periodic in the tangential direction which requires a periodic scalar tridiagonal 
solver for the system of finite difference equations in the computational plane. A Thomas 
tridiagonal solver is used for the non-periodic, normal direction. 
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For a given iteration, the model subroutine will require the input of Re, velocity 
components(u, v), and the physical grid geometry (x,y) from the flow solver. Using these 
informations, the subroutine will model the corresponding eddy viscosity of the flow. A 
flowchart of the Baldwin-Barth one-equation model subroutine is presented in Figure 6.1. 

6.2 Verification of the Baldwin-Barth Turbulence Model 

Testing was conducted to verify that the subroutine works properly. The constraints 
imposed by the "O" grid and conformal mapping complicated the verification process. The 
ideal test case would be a turbulent flow on a flat plate. With no pressure gradient, the flow 
characteristic is dominated by eddy viscosity since the molecular viscosity is small. How 
on a flat plate has been thoroughly investigated and analytical descriptions of the flow are 
available for comparisons. However, numerical computation of flow on a flat plate with 
"O" grid has proven to be extremely difficult The grid spacing of a zero-thickness, flat 
plate collapses near the surface and especially at the leading and trailing edges. Difficulties 
in calculating the transformation metrics makes proper transformation of such grids nearly 
impossible. A pseudo- flat plate was not used due to similar restrictions associated with 
conformal mapping. Therefore, a NACA-0012 airfoil was chosen as the test body. 

ARC2D and ZETA, using the Baldwin-Barth turbulence model, computed cases for 
comparison. The test cases ran at Re = 3.9 lx 10 6 and a = 5.0 degrees on a 120x80 grid. 
The velocity magnitude contours and eddy viscosity distributions of ZETA and ARC2D are 
shown in Figure 6.2 through Figure 6.5. Although the figures show global similar 
features, ZETA with the Baldwin-Barth model predicts flow separation early. Separation is 
indicated by flow reversal. As in the ARC2D case, ZETA with Baldwin-Lomax does not 
show any flow separation (Fig. 6.6). The corresponding vorticity is presented in Figure 
6.7. The unexpected flow separation could be caused by the incorrect implementation of the 
Baldwin-Barth model or undesirable interactions between ZETA and constraints imposed 
by the turbulence subroutine. The ZETA code and these probable causes were examined. 
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6.2.1 Test of the Isolated T urbulence Subroutine 

The field variable of the model subroutine is the turbulent Reynolds number, Rj 
which is directly proportional to eddy viscosity. Thus, the prediction of Rj must be correct. 
An effective way of testing the turbulence subroutine is to isolate and place it in a small test 
program. Receiving a converged ARC2D velocity solution field(u,v) along with the grid 
geometry(x.y), the subroutine would restart and iterate until Rj converged. 

The eddy viscosity, Vt, was calculated with the converged Rt of the test program 
and compared with the v t of ARC2D. Figure 6.8 shows the global features of the test case 
which are in good agreement with that of ARC2D (Fig. 6.5). The trailing edge features of 
the Vt distribution for both ARC2D and test case also compare well (Fig. 6.9 and Fig. 

6.10). A more detailed comparison of the Rj profiles at six locations on the more critical 
top surface of the airfoil is in Figure 6.1 1. The matching Vt distributions and Ry profiles 

indicate that the subroutine in ZETA predicts eddy viscosity correctly. 

Another convincing comparison is the velocity distributions u + vs y+, which is 
computed in the subroutine with the mean velocity and geometry, at mid-chord for both 
ZETA and ARC2D cases (Fig. 6.12). The good agreement of the velocity profiles at mid- 
chord as well as other chord stations also affirms that the model subroutine in ZETA is 
producing the correct eddy viscosity for a given flow. It is also interesting to note that the 
velocity profile is similar to that of a flat plate (Fig. 6.13). Although ZETA with the 
Baldwin-Barth model predicts flow separation early, the above tests show that the model 
subroutine is implemented correctly. The focus will be shifted to the ZETA code and its 
numerical interactions with the Baldwin-Barth model subroutine. 

However, it is important to note that ARC2D and the isolated model test cases are 
time-average runs but ZETA is time-accurate at all time. Recommended modifications for 
the time-accurate code were made to the model subroutine in ZETA. Since no known time- 
accurate computation with the Baldwin-Barth model is available, the model subroutine's 
numerical capability to compute time-accurately may require further examination. 


6.3 
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Omission of S© 

The governing equations of ZETA are the continuity equation and the vorticity 
transport equation. The vorticity transport equation (Equation 2.2) has neglected a scalar 
term, S w : 



This is justified below by the order-of-magnitude analysis of S^. 

For a laminar boundary layer, the boundary layer thickness is proportional to the 
square root of the kinematic viscosity [17]: 

8~Vv 

This boundary layer thickness is very small compared to the chord of the airfoil. The 
dimensionless 5 is therefore very small compared to unity. v e is of the order 5 as indicated 
in the numerical computations using both Baldwin-Lomax and Baldwin-Barth turbulence 
models. This implies that the inverse of the turbulent Reynolds number: 

X = -^L- 

R V„C 

is also of the order 8. Therefore, the orders of magnitude of the individual terms of S w , 
which is nondimensionalized with the free-stream velocity, Voo, and the chord, C, are as 
follow: 

S -J 32 f l9u ) d 2 ( ld v li d 2 ( ld y \ d 2 fld u | 

^ 3x 2 lRdy; dy2\Rdx) dxdy\Rdyl dxdy\Rdx) 

1 5 1 _L 5 5 1 5^ l 8 1 

8 8 2 8 8 8 

Thus, Soj is of order one. The vorticity transport equation is of the order 1/8. Therefore, 

the omission of S a in the vorticity transport equation is justified. 

Having shown that the turbulence model subroutine works properly and the 
derivation of the governing equations was correct, the numerical aspects of ZETA and its 
interactions with the turbulence model subroutine will then be investigated. 


38 

6.4 The Effects of Grid 

The grid has the utmost effect on computational accuracy of any case. Depending 
on the numerical methods of a code, constraints or demands on the grid will be different 
and sometimes conflicting. Grid requirements will also change with flow conditions. For 
conventional finite difference codes, a finer grid will generally increase the accuracy of the 
computation. This grid effect applies to the Baldwin-Barth model subroutine which utilizes 
the finite difference formulation. However, the integro-differential approach with Fourier 
series expansions of ZETA places unique demands on grid resolution. 

The Baldwin-Barth numerical formulation requires good grid resolution in both 
tangential and normal directions of the airfoil. To capture the flow features in the boundary 
layer requires proper grid spacing, especially in the flow direction and in the near-wall 
region. Having the appropriate first grid spacing off the airfoil surface is also essential. The 
Baldwin-Barth subroutine requires a grid point to be inside the laminar sublayer for proper 
estimation of the velocity profile [13]. The slope of the profile is used in computing the 
friction velocity in the subroutine. Inaccurate approximation of the slope will cause the 
calculation of y + , the damping factors, and hence eddy viscosity to be incorrect (Fig. 6.1). 
Computations with a grid having the first velocity grid point off the airfoil surface outside 
of the laminar sublayer will give inaccurate flow solutions. 

In contrast, the flow solver, ZETA, prefers a coarse grid with a minimum of only 
about 10 to 20 grid points inside the boundary layer [10]. Testings also found that the first 
radial grid point off the airfoil surface should be outside of the laminar sublayer. This 
unusual grid requirement is a result of the numerical formulation in the code. When the 
airfoil is impulsively started to move in the fluid, the kinematic aspect of the code computes 
boundary vorticity on the first ring of vorticity grid points off the surface (Fig. 2.2). This 
procedure integrates over the viscous zone outside of the first ring to obtain the boundaiy 
vorticity. Using the known velocities at the previous time level and the boundary 
vorticities, the vorticity transport equation computes the vorticity in the remaining viscous 
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zone. The boundary condition is updated at each time level. However, the kinematic aspect 
of the code could not compute constant vorticity at grid points inside the laminar sublayer. 
The numerical error will then accumulate as the time level advances and will have a great 
impact on the accuracy of the boundary vorticity and the overall flow solution. 

This premise is tested by running ZETA in the laminar mode where the effect of the 
turbulence model is removed. Changing the first radial velocity grid spacings off the 
surface has significant effects on the solution of the given flow (Fig. 6.14 to 6.16). The 
flow separation region progressively increases as the first grid spacing decreases from 
0.001 to 0.00006. The amount of separation should not change with the first radial grid 
spacing. The case with more grid points inside the laminar sublayer show unreasonably 
large amount of separations (Fig. 6.16). The numerical problem of updating the vorticity 
boundary inside the laminar sublayer appeared to cause the flow differences. Another 
probable reason for the separation differences is the higher numerical dissipation rate of 
coarse grids. Numerically, the flow is more likely to remain attached on coarser grid on 
which disturbances quickly dissipate. 

For the same flow condition and a first grid spacing of 0.0001, systematically 
reducing all the numerical eddy viscosity to 20%, 10%, and zero percent of the original 
Baldwin-Barth subroutine value does not change the flow features significantly (Fig. 6.17 
to 6. 19). The flows with reduced eddy viscosity behave like the normal turbulent flow. 
Such behavior suggests that the flow solver may be introducing some interferences with 
greater effects than the computed eddy viscosity. 

TFT A with the Baldwin-Barth model running on a grid with the first ring of grid 
points outside of the laminar sublayer did not show separation (Fig. 6.20 and 6.21) but the 
predicted aerodynamic loads are low. This is expected since the model requires a grid point 
within the laminar sublayer. However, having any grid point in the laminar sublayer will 
lead to inaccurate flow solution with ZETA. 

The conflicting grid requirements of ZETA and the Baldwin-Barth model make 
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accurate flow computation difficult unless a solution to the problem is found. For the above 
case using the Baldwin-Barth model (Fig. 6.2), a systematic investigation of grid size and 
first grid spacing showed that a grid of 120x80 with a first grid spacing of 0.0001 gives the 
best but still incorrect solution. Adding grid lines in the tangential direction did not improve 
the flow solution. On the same 120x80 grid, ZETA with the Baldwin-Lomax model under- 
predict the aerodynamic loads such as Cl by about 20% although it showed correct global 
features(Fig. 6.6). A coarser grid is recommended when using the Baldwin-Lomax 
turbulence model. A typical grid size of 80x45 was used for previous investigations of 
airfoils. The first grid spacing varies slightly with Re. 

6.5 The Effects of Reynolds Number 

Most numerical codes have a operating range of Reynolds numbers for which the 
resulting flow solutions are accurate. Generalizing from previous work, ZETA is not 
expected to perform well for cases where Re is less than 2x10 s . No other known Reynolds 
number effect on the flow solution is reported. Codes with the Baldwin-Barth model 
generally run at a Re of 2X10 6 or greater. Further investigation of Reynolds number effect 
was performed since the Baldwin-Barth model introduces new Re and grid demands. By 
running the NACA-0012 airfoil at a = 5° through a range of Reynolds numbers, it was 
found that Re affects ZETA with both Baldwin-Lomax and Baldwin-Barth models. 

Re from lxlO 4 to 4x10 s were tested with ZETA to examine the tangential velocity 
and vorticity along the top surface of the airfoil. The tangential velocity and vorticity along 
the airfoil surface start to fluctuate between Re of one and two millions with the Baldwin- 
Lomax model (Fig. 6.22 and 6.23). The normal velocity and vorticity exhibit no such 
fluctuations. There is also no fluctuation on the lower airfoil surface. The fluctuations did 
not disappeared when more lines were added in the tangential direction of the grid to 
improve the velocity gradient calculation. Cases with the Baldwin-Barth model showed 
more severe fluctuations in the tangential surface velocity and vorticity. Fluctuations begin 



at a Re of about lxlO 5 and becomes proportionally worse as Re increases (Fig. 6.24 and 
Fig. 6.25). The Baldwin-Barth model equation is a transport equation which may be the 
cause of the more severe fluctuations. Transport equations tend to accentuate disturbances 
in a code. Fluctuations also existed in the laminar mode using a fine grid. There is no 
evidence of fluctuations when using a coarse grids even at high Reynolds number (Fig. 
6.21). Having fluctuations in the laminar mode and with both turbulence models hint that 
the flow solver presently may not be able to accommodate a fine grid at higher Re- 
However, the grids used for ZETA could be considered coarse for most finite difference 
Navier- Stokes codes. 

6.6 The Performance of Turbulence Models 

For most applications of turbulence models including Baldwin-Lomax, Navier- 
Stokes codes would normally over-predict Cl- ZETA with the Baldwin-Lomax turbulence 
model however has been under-predicting Cl and the predicted Cl is generally lower at 
higher Re- The lift slope in Figure 6.26 was computed at Re = 3.91X10 6 and a = 5 degrees 
on a 80x45 grid. The experimental result on the figure is from Reference 2. The lower Cl 
values are consistent with ZETA's performance but the cause of the difference is not 
known. The performance of ZETA with the zero-equation model is generally acceptable for 
flows where separation is minimum. Flows with large separation region are common in 
unsteady aerodynamic investigations where dynamic stall is often encountered. As 
expected, ZETA's performance for these flows and high angle of attack flows are not as 
reliable. Two of the possible causes are compressibility effect and poor performance of the 
turbulence model in separated flows. The calculation of the length scale in the zero-equation 
model is known to be inaccurate for separated flows and this is true in ZETA. 

Though the overall performance of ZETA with the Baldwin-Lomax model is 
reasonable when used properly, it remains questionable. As shown in Figure 6.17 to 6.19, 
changing the eddy viscosity magnitude did not change the Baldwin-Barth flow solutions 


which are slightly separated When the Baldwin-Lomax model is applied, the flow remains 
attached Beside the eddy viscosity magnitude, the only other factor involved is the eddy 
viscosity distribution since the Baldwin-Barth subroutine can only affects the flow solution 
by the computed eddy viscosity. The eddy viscosity distribution of the Baldwin-Lomax 
case is typically like that of Fig. 6.7. Note that the eddy viscosity, v t , is constant in large 

area near the airfoil. This is unrealistic. Examining the velocity distribution, the Baldwin- 
Lomax case is also different from that of ZETA and ARC2D with the Baldwin-Barth model 
(Fig. 6.27). The reason for the attached flow when using the Baldwin-Lomax model is 
unclear. ZETA with the Baldwin-Lomax model however tends to delay flow separation 
(Fig. 6.26). 

There are advantages in using the Baldwin-Lomax model in ZETA. It is not time or 
memory intensive. A coarse grid(80x45) can be used and tends to provide better flow 
solutions than finer grids. With Baldwin-Lomax, steady and dynamic flow cases can also 
be ran quickly and interactively on the CRAY Y-MP. 

The Baldwin-Barth one-equation turbulence model has been implemented in mainly 
finite difference codes. The solution of the model equation requires a relatively significant 
addition of memory and time due to its added complexity and grid resolution requirement. 
Although the Baldwin-Barth model has been proven to be reliable in many flow conditions 
over airfoils, it has not been performing well in ZETA. As discussed in the previous 
sections, disruptive fine grids and Re factors might of contributed to the poor performance 
of ZETA with the Baldwin-Barth model. Therefore, extensive quantitative comparison with 
experimental results was not conducted. 

For a flow at Re = 3.91x106 and a = 20 degrees, the velocity contour shapes with 
both turbulence models are comparable (Fig. 6.28 and 6.29) but the values are different. At 
this condition, the airfoil is stalled and the surface boundary layer effects are negligible. 
Flow solvers are not expected to do well in these stall cases. However, the Baldwin-Barth 
turbulence model seems to perform better when the disruptive grid effect at higher Re is 


small. There is no fluctuation in the flow solution. The eddy viscosity distribution of the 
Baldwin-Barth model seems to be more reasonable than that of the Baldwin-Lomax model 
for separated flows (Fig. 6.30 and 6.31). One reason for the unrealistic eddy viscosity 
distribution of Baldwin-Lomax (Fig. 6.31) is that the subroutine imposed a maximum eddy 
viscosity limit of 100 times the laminar viscosity. In the development process [10], the 
Baldwin-Lomax model performed better with this imposed upper limit on eddy viscosity. 

6.6. 1 Aerodynamic Loads Prediction 

In the prediction of aerodynamic loads: Cl, Cd> and Cm, ZETA with Baldwin- 
Lomax performs well with a coarse, 80x45 grid (Fig. 6.26). Even the suction peak of the 
Cp distribution matches well with the experimental data (Fig. 6.32). The computation of C p 
however deteriorates with finer grids like those used for the Baldwin-Barth model. In the 
same figure, the C p distribution computed with the Baldwin-Barth model indicates that it 
performs well especially near the leading edge but shows early flow separation. The C p 
distribution fluctuates on the upper surface of the trailing edge. The fluctuation problem is 
most obvious in the time convergence history of Cl (Fig. 6.33). ZETA with the Baldwin- 
Lomax model shows no fluctuation on a coarse grid(80x45). The Cl compares well but it 
is expectedly lower than the experimental value of 0.58. It is much lower when using a 
finer grid(120x80). In the Baldwin-Barth case (Fig. 6.33), the average Cl = 0.48 is low 
but the peak Cl = 0.53 is comparable to the experimental Cl = 0.58. Cl does not fluctuate 
on grids with the first grid point outside of the laminar sublayer (dy = 0.01) but this model 
is not expected to perform well on such grid. The average of the fluctuating Cl values is 
also comparable to the experimental data before stall occurs (Fig. 6.34). However, ZETA 
with the Baldwin-Barth model is presently ineffective. The model would improve ZETA's 
performance in separated flows when the conflicting grid requirement is solved. Without 
the high frequency fluctuations, the results with the Baldwin-Barth model are promising. 


CHAPTER 7 


CONCLUDING REMARKS 


7.1 Conclusions 

The Baldwin-Barth one-equation turbulence model has been implemented 
into ZETA. The isolated turbulence model subroutine works correctly. ZETA with the 
Baldwin-Barth turbulence model can not be used for practical flow applications due to their 
conflicting grid requirements. Baldwin-Barth requires a point inside the laminar sublayer 
but ZETA does not work well with such grid. The implementation process revealed some 
limitations of ZETA and the integro-differential scheme. Major conclusions of this 
investigation includes: 

(1) The first grid point off the solid body surface should be outside of the 
laminar sublayer to compute accurately with ZETA. 

(2) For airfoils, ZETA should be run at Re greater than 2X10 5 . 

(3) A fine grid and higher Re will cause vorticity and velocity fluctuations. 

(4) The grid generator used for ZETA is simple but rigid. The inability to cluster 
grid lines and restriction to the use of the "O" grid hampered the study. 

(5) ZETA with the Baldwin-Lomax turbulence model seems to give better 
result, but the computed eddy viscosity distribution remains questionable. 

7.2 Recommendations 

Future development and work with ZETA should consider the use of an effective 
grid generator that is independent of the code and has adaptive capability. The conflicting 
grid requirement may be resolved by refining the computation of vorticity in the laminar 
sublayer in ZETA. The code should be expanded to include the compressible flow case. 
Recent work on vortex trapping to counteract unwanted dissipations has proven to be 
promising and ZETA should have such a feature since a zonal procedure is used. 
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TABLE I: Summary of Selected Test Cases 


4 5 


Study 

Case# 

Reynolds 

Number, 

Angle of 
Attack, a 

Grid Size 

First Grid 
Spacing, dy 

Code and 
Comments 

1 

3.91x106 

5° 

120x80 

0.0001 

Z/BB 

2 

3.91xl0 6 

5° 

120x80 

0.00006 

ARC2D 

3 

3.91X10 6 

5° 

120x80 

0.0001 

Z/BL 

B 

3.91x106 

5° 

120x80 

0.00006 

Isolated BB with 
ARC2D inputs 

1 5 

3.91x106 

5° 

120x80 

0.001 

Laminar 

1 6 

3.91x106 

5° 

120x80 

0.0001 

Laminar 

1 7 

3.91x106 

5o 

120x80 

0.00006 

Laminar 

8 

3.91x106 

5° 

120x80 

0.0001 

Z/BL with 20% of 
the computed v t 

9 

3.91x106 

50 

120x80 

0.0001 

Z/BL with 10% of 
the computed Vt 

10 

3.91x106 

5° 

120x80 

0.0001 

Z/BL with 0% of 
the computed v t 

11 

3.91x106 

50 

120x80 

0.01 

Z/BB 

12 

lxlO 6 

5° 

120x80 

0.0001 

Z/BL 

13 

2 x 106 

5° 

120x80 

0.0001 

Z/BL 

14 

lxlO 4 

5° 

120x80 

0.0001 

Z/BB 

15 

5x105 

50 

120x80 

0.0001 

Z/BB 

16 

2 x 106 

50 

120x80 

0.0001 

Z/BB 

17 

3.91x106 

3° 

80x45 

0.001 

Z/BL 

18 

3.91x106 

5° 

80x45 

0.001 

Z/BL 

19 

3.91x106 

10° 

80x45 

0.001 

Z/BL 

20 

3.91x106 

13° 

80x45 

0.001 

Z/BL 


Z/BB: ZETA with Baldwin-Barth(BB) turbulence model. 
Z/BL: ZETA with Baldwin-Lomax(BL) turbulence model. 
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TABLE I, Continued 
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Study 

Case# 

Reynolds 
Number, Re 

Angle of 
Attack, a 

Grid Size 

First Grid 
Spacing, dy 

Code and 
Comments 

21 

3.91x106 

15° 

80x45 

0.001 

Z/BL 

22 

3.91x106 

16° 

80x45 

0.001 

Z/BL 

23 

3.91xl0 6 

20° 

80x45 

0.001 

Z/BL 

24 

lxlO 6 

200 

120x80 

0.0001 

Z/BB 

25 

1x106 

20° 

120x80 

0.0001 

Z/BL 


Z/BB: ZETA with Baldwin-Barth(BB) turbulence model. 
Z/BL: ZETA with Baldwin-Lomax(BL) turbulence model. 





































Figure 1.1: UH-60A Black Hawks in forward flight. 
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Figure 3.1: Clustering of radial lines at the leading and trailing edges. 
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Figure 6.1: Flowchart of the Baldwin-Barth subroutine. 













Figure 6.2: Velocity contours of an NACA-0012 airfoil at R e =3.91x10^ 
and a = 5° from ZETA with Baldwin-Barth turbulence model. 



Figure 6.3: Eddy viscosity contours of an NACA-0012 airfoil at Re =3.91xl0 6 
and a = 5° from ZETA with Baldwin-Barth turbulence model. 




54 



Figure 6.4: Velocity contours of an NACA-0012 airfoil at Re =3.91x106 

and a = 5° from ARC2D with Baldwin-Barth turbulence model. 



Figure 6.5: Eddy viscosity contours of an NACA-0012 airfoil at Re =3.91x106 
and a = 5° from ARC2D with Baldwin-Barth turbulence model. 





Figure 6.6: Velocity contours of an NACA-0012 airfoil at =3.91xl0 6 

and a = 5° from TFT A with Baldwin-Lomax turbulence model. 



Figure 6.7: Eddy viscosity contours of an NACA-0012 airfoil at =3.91x10^ 
and a = 5° from TF T A with Baldwin-Lomax turbulence model. 
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Figure 6.8: Eddy viscosity contours of an NACA-0012 airfoil. at Re =3.91x10^ 

and a = 5° from isolated Baldwin-Barth turbulence model test program. 



Figure 6.9 : Trailing edge eddy viscosity contours of an NACA-0012 airfoil at 

Re =3.9 lxl 0 6 and a = 5° from ARC2D with Baldwin-Barth model. 



Figure 6.10: Trailing edge eddy viscosity contours of an NACA-0012 airfoil at 

Re -3,91xl0 6 and a = 5° from the isolated Baldwin-Barth subroutine. 
























Figure 6. 14: Velocity contours with first grid spacing, dy = 0.001. 



Figure 6.15: Velocity contours with first grid spacing, dy = 0.0001. 







Figure 6.18: Velocity contours with 10% of computed eddy viscosity 
from the Baldwin-Barth turbulence model. 
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Figure 6. 19: Velocity contours with 0% of computed eddy viscosity 
from the Baldwin-Barth turbulence model. 
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Figure 6 . 20 : Velocity contours of an NACA-0O12 airfoil at Re =3.91x106 and 

a = 5° from ZETA with Baidwin-Barth turbulence mode, dy = 0.01 




Figure 6.21: Eddy viscosity contours of a NACAQ012 at R« *3.91x106 and 

a = 5° from ZETA with Baidwin-Barth turbulence model, dy = 0.01 














Figure 6.24: Surface velocity flucturations of ZETA with Baldwin-Barth 
turbulence model at different R^. dy = 0.0001 
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Figure 6.27: u + vs y + at mid-chord of ARC2D with Baldwin-Barth 
turbulence model and ZETA with Bladwin-Barth and 
Bladwin-Lomax turbulence models. 
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Figure 6.28: Velocity contours of an NACA-0012 airfoil at Re = lxlO 6 

and a = 20° from ZETA with Baldwin-Barth turbulence model. 
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Figure 6.29: Velocity contours of an NACA-0012 airfoil at Re = 1x10^ 

and a = 20° from ZETA with Baldwin-Lomax turbulence model. 
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Figure 6.30: Eddy viscosity contours of an NACA-0012 airfoil at Re =1x10^ 
and a = 20° from ZETA with Baldwin-Barth turbulence model. 
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Figure 6.31: Eddy viscosity contours of an NACA-0012 airfoil at Re =1x10^ 
and a = 20° from ZETA with Baldwin-Lorn ax turbulence model. 




Figure 6.32: Computed and experimental pressure coefficients of 
an NACA-0012 airfoil at R e =3.91x10^ and a = 5°. 








APPENDIX A 


LISTING OF THE BALDWIN-BARTH SUBROUTINE 


£********************************************** * *********************** 
C*************** turbulent viscosity ******************************** 
£****♦ ******************* ******** * *** ********************************** 
SUBROUTINE BBEDDY 

C ONE EQUATION TURBULENCE MODEL 
C VERSION 1.0 

C BY BALDWIN AND BARTH NASA AMES R.C. 

C ********* GENERAL note *********** 

C CONSTANTS ARE GIVEN BELOW FOR THE MODEL THAT APPEARED IN 

C NASA TM-102847. THE SECOND VERSION THAT APPEARED IN 
C AIAA PAPER 91-0610 IS NOT USED. 

PARAMETER(JDIM=120XDIM=80) 

PARAMETER(MAXJ=121,MAXK=81) 

COMMON /GRD/JMAXJM2XFC.KMAXN 

COMMON /BB_1/TURMU(JDIM,KDIM),TURRE(JDIM,KDIM),RE .FIRST 
COMMON /BB_2/ SMIN(MAXJ,MAXK), BWT(MAXJ,MAXK), 

$ IB(MAXJ,MAXK),X(JDIM,KDIM),Y(JDIM,KDIM), 

$ XY(JDIM,KDIM.4)XYJ(JDIM,KDIM). 

$ AA(0:MAXJ,0:MAXK),AAN(MAXJ,MAXK) 

COMMON/BBJ3/U(MAXJ,MAXK),V(MAXJ,MAXK)pSMACH,ALPI,T,AL 
COMMON/DELT A/DZ .DTET JDT 

COMMON /TUR l/USTAR(JDIM)£DDY(JDIM,KDIM)JOTUR(JDIM), 

$ YN(MAXJ JCDIM) JOT JQTB 

COMMON /COE/AF.UI.VI.OMG.VSCNPL.ICTUR.ICST.ICPL 

DIMENSION Q(JDIM,KDIM,4),VORT(JDIM,KDIM) 

DIMENSION PRESS(JDIM,KDIM) 

DIMENSION FNU(JDIMJCDIM)T)S(JDIMJCDIM) 

DIMENSION mXMAXJ,MAXK),VD(MAXJ.MAXK),QD(JDIM,KDIM), 
DIMENSION FND(KDIM),TEMPU(MAXJ,MAXK), TEMP V (M AXJ ,MAXK) , 

> EMPS(MAXJ,MAXK), TEMPMU(JDIM,KDIM), 

> GQ(0:JDIM,0:KDIM,2,2), 

> GQN(JDIM,KDIM,2,2), 

> DAMP 1 (JDIM.KDIM), DAMP 1 M(JDIM,KDIM), 

> DAMP2(JDIM,KDIM), WORKX(MAXKMAXJ), 

> WORRY (MAXJ .MAXK) 

DIMENSION AX(KDIMJDIM)3X(KDIMJDIM).CX(KDIMJDIM). 

> DX(KDIMJDIM)£X(KDIMJDIM)PX(KDIMJDIM). 

> AY(JDIM,KDIM)3Y(JDIM,KDIM) 1 CY(JDIM,KDIM). 

> DY(JDIM,KDIM)3Y(JDIMJCDIM)JY(JDIMJCDIM) 

LOGICAL FIRST 
DATA FIRST /.TRUE. / 


77 


78 

C THE FOLLOWING TWO VARIABLES ARE FOR 'O' GRID. 

C THEY ARE 2 AND JMAX-1 FOR 'C' GRID. 

JTAIL1 = 1 
JTAIL2 = JMAX 

C THESE ARE THE CONSTANTS FOR THE MODEL APPEARING IN NASA 
C TM- 102847. 

AKARMAN = .41 
CMU = .09 
C1E = 1.20 
C2E = 2 00 

SIGE = AKARMAN** 2/((C2E-C 1 E)* SQRT (CMU)) 

APLUS1 = 26. 

APLUS2 = 10. 

RETINF = IP-12 

C DTM IN THE MODEL IS EQUAL TO DT OF THE SOLVER AND NNIT IS 
C ONE FOR TIME ACCURATE SOLUTION. 

NNIT = 1 
DTM = DT 

C FOR STEADY-STATE CALCULATIONS WE USUALLY TAKE A TIME STEP 
C IN THE TURBULENCE MODEL WHICH IS LARGER THAN THE FLOW 

C SOLVER (BUT NOT TOO LARGE). 

C IF(ICST .EQ. 0 .OR. ICST .EQ. 1)DTM = 200 *DT 

C IF(DTM.GT.50)DTM = 50. 

C **THESE CONSTANTS FOR THE SUTHERLAND'S LAW ARE NOT USED 
C IN THE INCOMPRESSIBLE FLOW 

C C2B=198.6/TINF 

C C2BP = C2B + 1. 

IF(FTRST)THEN 

FIRST = PALSE. 

NNIT =1 

WRITE(6 ,*)' KARMAN CONSTANT = '.AKARMAN 
WRITE(6 *)' CMU = '.CMU 

WRITE(6 *)' C1E = '.C1E 

WRITE(6 *)' C2E = '.C2E 

WRITE(6 *)' SIGMA E = '.SIGE 

WRITE(6,*)' APLUS 1 = \APLUS 1 

WRITE(6,*)' APLUS2 = ’.APLUS2 

WRITE(6,*)' RE_T INF = ’.RETINF 

WRITE(6,*)' RE = 'PE 

WRITE(6 ,*)' JDIM.KDIM = ’ JDIM.KDIM 

WRITE(6 ,*)' JMAX, KM AX = 'JMAX.KMAX 

WRITER,*)' MAXJ.MAXK = '.MAXJ.MAXK 

WRITE(6 ,*)' JTAIL1 JTAIL2 = 'JTAIL1 JTAIL2 

WRITE(6,*)' DTM = '.DT 

IF(JTAIL2 .EQ. JMAX)WRITE(6,*)' JTAIL2 = JMAX ' 



c 


COMPUTE METRICS 
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CALL XYMETS(JMAX,KMAX,X,Y XY ,XYJ) 

DO J= UMAX 
DO K=1,KMAX 

C ****FOR TIME ACCURATE SOLUTION**** 

DS(JJC) = 1. 

C ****FOR NON-TIME ACCURATE SOLUTION**** 

C DS(JJC) = l./(l. + SQRT(XYJ(J,K))) 

ENDDO 

ENDDO 

C COMPUTE GENERALIZED DISTANCE FUNCTION WHICH IS 
C MINIMUM DISTANCE TO WALL. 

C SAME DISTANCE FUNCTION APPEARS IN AIAA PAPER 91-0721 


DO J= UMAX 
DO K=1,KMAX 

DISTMS = 10000. 

DISTML = 10000. 

SMIN(JJC) = 10000. 

IB(JJC) =0 

DO JB=JTAIL1 JTAIL2-1 
JJ1 = JB 
JJ2 = JB + 1 

DXA = X(JJ2,1) - X(JJ1,1) 

DYA = Y(JJ2,1) - Y(JJ1,1) 

SNX = -DYA 

SNY = DXA 

DXX1 = X(JJC>X(JB,1) 

DYY1 = Y(JJC)-Y(JB,1) 

DXX2 = X(JJC)-X(JB+1,1) 

DYY2 = Y(J JC)- Y (JB+ 1,1) 

DXX = .5*(DXX1+DXX2) 

DYY = ,5*(DYY1+DYY2) 

DIST1 = (DXX1**2+DYY1**2) 

DIST2 = (DXX2* *2+DY Y2* *2) 

DISTS = MIN(DISTU)IST2) 

DISTL = MAX(DIST1 .DIST2) 

DOT = SNX*DXX + SNY* DYY 

IF((DISTSXE.DISTMS).AND.(DOT.GT.-1.E-10))THEN 

DISTMS = DISTS 

DISTML = DISTL 

SS 

> = (DXA*(X(J JC)-X(JJ 1,1)) 

> + DYA*(Y(JJC)-Y(JJ1,1)))/(DXA**2+DYA**2) 

SS = M AX(0.0,MIN(S S , 1 .0)) 

XPT = X(JB,1) + SS*(X(JB+1,1)-X(JB,1)) 

YPT = Y(JB,1) + SS*(Y(JB+1,1)-Y(JB,1)) 

SSMIN = SQRT((X(J JC)-XPT)**2+(Y (J,K)-YPT)**2) 
IF(SSMIN XT. SMIN(JX))THEN 
SMIN(J,K) = SSMIN 
BWT(J,K) = SS 
IB(JX) = JB 
DISTMS = DISTS 
DISTML = DISTL 
ENDIF 
ENDIF 



nnnnnonn 
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ENDDO 

ENDDO 

ENDDO 

C COMPUTE WEIGHING FACTORS FOR VELOCITY GRADIENT 

C CALCULATION. 


DO K=1.KMAX 
DO J= UMAX 
AA(JJC) = 0.0 
AAN(J.K) = 0.0 
ENDDO 
ENDDO 


DO J= UMAX 
DO K=1 JCMAX-1 
DXA = X(J,K+1)-X(JJC) 

DYA = Y(J,K+1)-Y(JJC) 

SNX = DYA 
SNY = -DXA 
JM1 = J-l 

IF (J .EQ. 1) JM1=JMAX 

AA(JM1,K) = AA(JM1,K) + SNX* ,5*(X(J,K+ 1 )+X(J,K)) 
AA(JJC) = AA(JJC) - SNX*.5*(X(J,K+1>+X(J,K)) 
ENDDO 
ENDDO 


DO J= UMAX 
JP1=J+1 

IF(J .EQ. JMAX) JP1 = 1 
DO K=1,KMAX 
DXA = X(JPUC)-X(J,K) 

DYA = Y(JP1 JC)-Y(J,K) 

SNX = DYA 
SNY = -DXA 

AA(JX-l) = AA(JX-1) - SNX*.5*(X(JP1,K)+X(UC)) 
AA(J,K) = AA(J,K) + SNX* .5* (X(JP 1 ,K)+X(J,K)) 
ENDDO 
ENDDO 
DO J= UMAX 
JP1=J+1 

IF(J JEQ. JMAX) JP1 = 1 
DO K=1 JCMAX-1 
AAN(JJC) = AAN(J,K) + AA(JJK) 

AAN(JP1 ,K) = AAN(JP1,K) + AA(J,K) 

AAN(J,K+1) = AAN(J,K+1) + AA(J,K) 
AAN(JPUC+1) = AAN(JP1,K+1) + AA(JJC) 

ENDDO 

ENDDO 

IF(JTAIL2 .NE. JMAX)THEN 
K = 1 

DOJ-UTAIL1 
JJ = JMAX - J + 1 
TEMPI = AAN(J.K) 

TEMP2 = AAN(JJ,K) 

AAN(J,K) = TEMPI + TEMP2 
AAN(JJJC)= TEMP2 + TEMPI 



on 
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ENDDO 
ENDIF 

IF(ICST .EQ. 0)THEN 

PRINT* ,’********INmAL RUN********' 

DO K=1,KMAX 
DO J= UMAX 
TURRE(JJC) = RETINF 
ENDDO 
ENDDO 
ELSE 
NNIT = 1 

PRINT* ,'*******REST ART********' 

OPEN(16) 

READ(16 *) JAJCA 

READ(16 ,*) ((TURRE(JJC), J=UA), K=1,KA) 

REWIND 16 
ENDIF 

WRITE(* ,*)' FINISHED TURB MODEL INIT ' 

ENDIF 

C DO J=JTAIL1 JTAIL2 

C TURRE(J,1) = 0.0 

C ENDDO 


JX=JMAX+1 
DO K=1,KMAX 
SMIN(JX.K)=SMIN(1 X) 
ENDDO 


C ********CHANGE J DIRECTION FOR VELOCITIES********* 

DO K = 1, KMAX 
DOJ = 2, JMAX 
TEMPU(JX) = U(JX) 

TEMPV(JJC) = V(JJC) 

ENDDO 

ENDDO 

C 

DO K = 1, KMAX 
DOJ = 2, JMAX 
UD(JJC) = TEMPU(JMAX+2-J,K) 

VD(J,K) = TEMPV (JMAX+2-J JC) 

ENDDO 

UD(UC) = U(UC) 

VD(UC) = V(UC) 

ENDDO 

C 

C NOW THE MODEL 

C COMPUTE THE VELOCITY GRADIENT VECTOR 


DO K= 1 ,KMAX 
DO J=1JMAX 
CKJJC,1)=1. 
Q(J,K,2) = UD(J,K) 
Q(J,K,3) = VD(JX) 
FNU(J,K) = 1. 
ENDDO 
ENDDO 
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C THIS ROUTINE ACTUALLY COMPUTES THE GRADIENT OF THE MEAN 

C VELOCITY 

CALL GR AD V(M AXJ ,M AXK JDIM.KDIM JMAX ,KMAX JT AIL 1 JTAIL2 , 

> X, Y ,A A, AAN,Q,GQ,GQN) 

C GENERAL NOTE: THE EDDY VISCOSITY IS SUPPLIED AT HALF POINTS 

C IN K, NU(J,K+ 1/2). THIS MEANS THAT SOME QUANTITIES ARE NEEDED 

C AT HALF POINTS AND OTHERS ARE NOT. 

C STEP 1. CALCULATE THE DAMPING FACTORS AT 1/2 POINTS 

RATT= C1E/C2E 

DO J=1 JMAX 
DO K=1,KMAX-1 
JB = EB(J,K) 

WT = BWT(JJC) 

JJ1 = JB 
JJ2 = JB + 1 

IF(JB .EQ. JMAX) JJ2 = 1 
WM1 = FNU(JJ1,1) 

WM2 = FNU(JJ2,1) 

SSW1 = ABS(GQ(JJ1,U,2) - GQ(JJ1,1,2,1)) 

SSW2 = ABS(GQ(JJ2, 1,1,2) - GQ(JJ2, 1,2,1)) 

SSW = ABS(SSW1 + WT*(SSW2-SS W 1)) 

WNU = WM1 + WT*(WM2-WM1) 

RHOW = Q(JJ1,1,1) + WT*(Q(JJ2,1,1)-Q(JJ1,1,1)) 

RA = SQRT( RE*(SSW/WNU + 1 .E-10*FSMACH**2)) 

YPLS = RA*SMIN(J,K) 

YPLS = MAX(YPLS,.0001) 

YPLSM = RA*.5*(SMIN(J,K)+SMIN(J,K+1)) 

IF(J .EQ. 91 .AND. K .LE. 5) PRINT*, 'Y+= '.YPLSM 
EXP1 = EXP(-(YPLS/APLUS 1)) 

EXP2 = EXP(-(YPLS/APLUS2)) 

DAMP1(J,K) = (l.-EXPl)*(l.-EXP2) 

DAMPIM(JJC) = (1. - EXP(-(YPLSM/APLUS 1)))* 

> (1 . - EXP(-(YPLS M/APLUS2))) 

DDY = EXP1/APLUS 1*(1 .-EXP2) 

> + EXP2/APLUS2*(1.-EXP1) 

DD = DAMP1(J,K) 

SDD = SQRT(DD) 

DAMP2(J,K) = 

> RATT + (1 .-RATT)*( 1 ,/(AKARMAN*YPLS) + DD )* 

> (SDD + YPLS *DDY/SDD) 


ENDDO 

DAMP1 (JJCMAX) = 1.0 
DAMP 1 M(J,KMAX) = 1.0 
ENDDO 




NOW SOLVE THE EQUATION 
DO 500NIT=1 J'TNIT 


DO K=1,KMAX 
DO J= 1JMAX 

TURMU(UC) = CMU*DAMP1 (J,K)*TURRE(J,K) 
ENDDO 
ENDDO 

F_ETA_ETA VISCOUS TERMS 


DO J=1 JMAX 
DO K=2,KMAX-1 
KP1 = K+l 
KM1 =K-1 

XY3P= .5*(XY(J,K3)+XY(J,KP1,3)) 

XY4P= .5*(XY(J,K,4)+XY(J,KP1,4)) 

TIP = (XY3P*XY(J JC,3)+XY4P*XY(J,K,4)) 

XY3M = .5*(XY(J,K3)+XY(J ( KM1,3)) 

XY4M = ,5*(XY(J,K,4)+XY(J,KM1,4)) 

TTM = (XY3M*XY(J,K,3)+XY4M*XY(J,K,4)) 

CNUD = ( FNU(JJC)+TURMU(J,K)/SIGE 

+ CMU*DAMP1(J,K)*TURRE(JJC)/SIGE )/RE 

CDP = TTP*CNUD 
CDM = TTM*CNUD 

TREP = ,5*(TURRE(J,KP1)+TURRE(J,K)) 

TREM = ,5*(TURRE(UCM1)+TURRE(J,K)) 

CAP = CMU*TREP*DAMP1M(J,K )*TIP/(SIGE 1, 'RE) 
CAM = CMU*TREM* DAMP 1 M(J JCM 1 )*TTM/(SIGE*RE) 

THIS COMES FROM MAXIMUM PRINCIPLE ANALYSIS 

BY(J,K) = MIN(-CDM + CAM.0.0) 

CY(UC) = MAX( CDP - CAP.0.0) + MAX( CDM - CAM.0.0) 
DY(UC) = MIN(-CDP + CAP.0.0) 

FY(JJC) = - BY(JJC)*TURRE(J,KMI) 

- CY(J,K)*TURRE(J,K ) 

- DY(J,K)*TURRE(J,KP1) 


ENDDO 

ENDDO 

ADVECTIVE TERMS IN ETA 


DO J= UMAX 
DO K=2,KMAX-1 

UU = (XY(UC,3)*Q(J,K.2)+XY(J > K,4)*Q(J,K.3)) 
SGNU = SIGN(1.,UU) 

APP = .5*(1.+SGNU) 

APM = .5*(1.-SGNU) 

FY(J,K) = FY(UC) - 


UU*( APP* (TURRE(J JC )-TURRE(J,K- 1 )) 
+APM*(TURRE(J,K+1)-TURRE(J,K)) ) 
BY(J,K) = BY(J,K) - UU*APP 
CY(J,K) = CY(J,K) + UU* ( APP- APM) 

DY(J,K) = DY(J,K) + UU*APM 
ENDDO 
ENDDO 

PRINT*, "ETA TERMS DONE' 

E_XI_XI VISCOUS TERMS 


DO J= UMAX 
JP1 = J+l 
JM1 = J-l 

EF(J TQ. 1) JM1 = JMAX 
IF(J .EQ. JMAX) JP1 = 1 
DO K=2,KMAX-1 

XY1P = ,5*(XY(J,K,1)+XY(JP1,K,1)) 

XY2P= .5*(XY(J,K,2)+XY(JP1,K,2)) 

TTP = (XY1P*XY(JJC,1)+XY2P*XY(UC,2)) 

XY1M= ,5*(XY (J,K, 1 )+X Y (JM1 ,K, 1 )) 

XY2M = .5*(XY(J,K,2)+XY(JM1,K,2)) 

TTM = (XY 1 M*XY(J,K, 1 )+XY2M*XY(J,K,2)) 

CNUD=( FNU(J,K)+TURMU(J JC)/SIGE 

+ CMU*DAMP1(JJC)*TURRE(UC)/SIGE )/RE 


CDP = TTP*CNUD 
CDM = TTM*CNUD 

TREP =.5*(TURRE(JP1 ,K)+TURRE(JJC)) 

TREM =.5*(TURRE(JMUC)+TURRE(JJ0) 

CAP=CMU*TREP*.5*(DAMP1(J,K)+DAMP1(JP1,K))*TTP/(SIGE*RE) 
CAM=CMU*TREM* .5*(DAMP 1 (JM 1 ,K)+DAMP1(J,K))*TTM/(SIGE*RE) 


TfflS COMES FROM MAXIMUM PRINCIPLE ANALYSIS 
BX(K J) = MIN(-CDM + CAM.0.0) 

CX(KJ) = MAX( CDP - CAP.0.0) + MAX( CDM - CAM.0.0) 
DX(KJ) = MIN(-CDP + CAP, 0.0) 

FY(J,K) = FY(J,K) - BX(KJ)*TURRE(JM1,K) 

- CX(K J)*TURRE(JJC ) 

- DX(K r J)*TURRE(JP 1 ,K) 


ENDDO 

ENDDO 

ADVECTIVE TERMS IN XI 


DO J= UMAX 
PI = J+l 
JM1 = J-l 

IF(J JEQ. 1) JM1 = JMAX 
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IF(J .EQ. JMAX) JP1 = 1 
DO K=2,KMAX-1 

UU = (XY(JJC,1)*Q(J,K^)+XY(J,K ( 2)*Q(J,K,3)) 
SGNU = SIGN(1.,UU) 

APP = .5*(1.+SGNU) 

APM = ,5*(1.-SGNU) 

FY(J,K)= FY(J,K) - 

> UU*(APP* (TURRE(J ,K>TURRE(JM 1 ,K)) 

> +APM*(TURRE(JP1,K)-TURRE(JJC)) ) 
BX(KJ) = BX(KJ) - UU*APP 

CX(KJ) = CX(KJ) + UU*(APP-APM) 

DX(K J) = DX(KJ) + UU*APM 
ENDDO 
ENDDO 

PRINT*, XI TERMS DONE’ 

C 

DO K=2,KMAX-1 
KM1 = K - 1 
DO J=1 JMAX 
UX = GQN(JJC,1,1) 

UY = GQN(JJC,U) 

VX = GQN(JJC,2,1) 

VY = GQN(JJC,2^) 

C REAL PRODUCTION OF K 

C SS = SQRT(ABS( 

C > 2*(UX**2+VX*UY+VY**2)+UY**2+VX**2 

C > -.666666*(UX + VY)**2)) 


C THIN LAYER APPROXIMATION 


SS = ABS(UY-VX) 


TT=(C2E*DAMP2(J,K>C1E)*SQRT(CMU*DAMP1(JJC))*SS 

FACT = DS(JJC)*DTM 
BX(K J) = BX(KJ)*FACT 
CX(KJ) = CX(KJ)*FACT + 1. 

DX(K J) = DX(KJ)*FACT 
BY(JX) = BY(JJC)*FACT 
CY(JJC) = CY(JJC)*FACT + 1. 

DY(JX) = DY(J,K)*FACT 
FY(J,K) = (FY(J JC)+TT*TURRE(JJC))*FACT 
ENDDO 
ENDDO 

CALL TRTV(JDIM,KDIM,1 JMAX.2.KMAX- 1 ,B Y.CY.DY JFY) 


DO K=2,KMAX-1 
DO J= UMAX 
BY(UQ = BX(KJ) 
CY(JJC) = CX(KJ) 
DY(J,K) = DX(KJ) 
ENDDO 
ENDDO 


nnnonn 


> 


CALL VTRIP(JDIM,KDIM, 1 JMAX,2,KMAX- 1 , 
BY.CY.DYJT) 
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DO K=2,KMAX-1 
DO J= UMAX 
FX(KJ) = FY(J,K) 
ENDDO 
ENDDO 
NEGN = 0 
SUMN = 0.0 
JMIN = 1 
KMIN = 1 

FXM = ABS(FX(1,1)) 


DO K=2,KMAX-1 
DO J= 1JMAX 

TURRE(J,K) = TURRE(J,K) + FX(K J) 
IF(ABS(FX(KJ)).GT.FXM)THEN 
FXM = ABS(FX(K J)) 

MIN = J 
KMIN = K 
ENDIF 

C TfflS NEXT CHECK IS NEVER NEEDED IN PRACTICE 
IF(TURRE(JJC)1T. l.E-12)THEN 
NEGN = NEGN + 1 
TURRE(J,K) = l.E-12 
FX(K J) = 0.0 
ENDIF 

SUMN = SUMN + (FX(KJ))**2 
ENDDO 
ENDDO 

SUMN = SQRT(SUMN)/FLOAT((JMAX-l)*(KMAX-l)) 

OUT FLOW BOUNDARY CONDITION FOR A POLAR MESH. 

DO J= 1JMAX 
TURRE(J.l) = 0. 

IF(Q(J JCMAX,2)*X Y (J JCMAX,3)+Q(J,KMAX,3)*XY(J,KMAX,4) .GE. 0.) 
* TURRE(J,KMAX) = TURRE(J,KMAX-1) 

ENDDO 

XMU=0.0 
TMAX = 0.0 
DO K=1,KMAX-1 
DO J= UMAX 

TMAX = MAX(TMAX,TURRE(J X)) 

TURMU(JJC) = CMU*DAMP1M(J,K)*.5*( 

> TURRE(J,K)*Q(UC,1) 

> + TURRE(J,K+1)*Q(UC+1,1)) 

IF (TURMU(JX) .GT. XMU) THEN 
XMU = TURMU(J,K) 

JMU=J 
KMU=K 
ENDIF 
ENDDO 
ENDDO 
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WRrrE(*,’(A^12.6,A^12.6^£12.6J5,lX45)') 

$ ' RESID RE_T = \SUMN,' MAX RT = \TMAX, 

$ ’RESID MAX = , ,FXM r JMIN,KMIN 

WRITER,*) ’MAX TURMU = XMU,XMU*AL/RE,’ (J,K)= ’JMU.KMU 

500 CONTINUE 

c ********C H ANGE J DIRECTION OF TURMU AND SET TO EDDY********* 

DO K = 1, KM AX 
DO J = 2, JMAX 
TEMPMU(J.K) = TURMU(J.K) 

ENDDO 

ENDDO 

C 

DO K = l.KMAX 
DO J = 2, JMAX 

EDDY(J,K) = TEMPMU(JMAX+2-J,K)*AL/RE 
ENDDO 

EDDY(1 ,K) = TURMU(1,K)*AL/RE 
ENDDO 
DO K= l.KMAX 

EDDY(1JC)=.5*(EDDY (2JQ+EDDY (JMAXJC)) 

ENDDO 


DO J = 1, JMAX 
EDDY(J,1) = 0.0 
ENDDO 

RETURN 

END 


£********************************************************************* 
C*******************^^ VELOCITY GRADENT*********************** 

C********************************************,', ** * ***** * * *** *** * + ** * ** + 

SUBROUTINE GRADV(MAXJ,MAXKJDIM,KDIMJMAX,KMAX, 

> JTAIL1 ,JTAIL2,X,Y,AA,AAN,Q,GQ,GQN) 

DIMENSION AA(0:MAXJ,0:MAXK),AAN(MAXJ,MAXK) 

DIMENSION Q(JDIM,KDIM,4) 

DIMENSION X(JDIM,KDIM),Y (JDIM.KDIM) 

DIMENSION GQ(0:JDIM,0:KDIM, 2, 2), GQN(JDIM,KDIM,2,2) 

C NOTE THE ZERO STARTING INDEX IN THE DIMENSIONS 

DO10N=U 
DO 10 K=1 JCMAX 
DO 10 J= 1JMAX 
GQ(J,K,N,1) = 0.0 
GQ(J.K,N,2) = 0.0 
GQN(J,K,N,1) = 0.0 
GQN(J,K,N,2) = 0.0 
10 CONTINUE 


DO 20 NN=2,3 
N = NN-1 

DO 20 K=1JCMAX-1 
DO 20 J= UMAX 
JM1 = J-l 

E(J EQ. 1) JM1 = JMAX 
DX = X(J,K+1)-X(JJK) 

DY = Y(J,K+1)-Y(J,K) 

SNX = DY 
SNY = -DX 

GQ(JM1,K,N,1) = GQ(JM1 + SNX*.5*(Q(J,K+1,NN)+Q(J,K,NN)) 
GQ(J,K,N,1) =GQ(JJC^,1) - SNX*.5*(Q(J,K+1,NN)+Q(J,K,NN)) 
GQ(JM1,K,N,2) = GQ0Ml < KJf l 2) + SNY* ,5*(Q(J,K+ 1 ,NN)+Q(J,K,NN)) 
GQCJJCJsT^) = GQ(J,K,N,2) - SNY*.5*(Q(J,K+1,NN)+Q(J > K,NN)) 

20 CONTINUE 


DO 30 NN=2,3 
N = NN-1 

DO 30 K= 1 ,KMAX 
DO 30 J= UMAX 
P1=J+1 

E(J .EQ. JMAX) El = 1 
DX = X(JP1 JC)-X(JJC) 

DY = Y(JP1,K)-Y(J,K) 

SNX = DY 
SNY = -DX 

GQ(J,K-1,N,1) = GQ(J,K-1,N,1) - SNX*.5*(Q(E1,K^1N)+Q(J,KJn1N)) 
GQ(J,K,N,1) = GQ(J,K,N,1) + SNX*.5*(Q(E1,K,NN)+Q(J,K,NN)) 
GQ(J,K-1^) = GQ(J,K-1,N,2) - SNY* .5* (Q(JP 1 ,K,NN)+Q(J,K,NN)) 
GQ(JJC,N,2) = GQ<J,K,N,2) + SNY*.5*(Q{EUCJ4N)+Q(J,KJ4N)) 
30 CONTINUE 
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DO 90 M=l,2 
DO90N=U 
DO 90 K=1,KMAX-1 
DO 90 J= UMAX 
JP1=J+1 

IF(J .EQ. JMAX) JP1 = 1 
GQN(J,K,N,M) = GQN(J,K,N,M) + GQ(J,KJM,M) 
GQN(JP1,K,N.M) = GQN(JP1,K ( N,M) + GQ(J,K,N,M) 
GQN(J.K+1,N > M) = GQN(J,K+1,N,M) + GQ(J,K,N,M) 
GQN(JP1,K+1,N,M) = GQN(JP 1 JC+ 1 ,N,M) + GQ(J,K,N,M) 
90 CONTINUE 

IF(JTAIL2 .NE. JMAX)THEN 
K= 1 

D0 95N=U 
DO 95 J=1 JTAIL1 
JJ = JMAX - J + 1 
TEMPI = GQN(J,K,N,1) 

EMP2 = GQN(JJ*K,N,1) 

GQN(J,K,N,1) = TEMPI + TEMP2 
GQN(JJJCJvT,l)= TEMP2 + TEMPI 
TEMPI = GQN(J,K,N,2) 

TEMP2 = GQN(JJJC,N^) 

GQN(J JC.N.2) = TEMPI + TEMP2 
GQN(JJ,K,N,2)= TEMP2 + TEMPI 
C95 CONTINUE 
C ENDIF 

DO110N=U 
DO 110 K=1,KMAX-1 
DO 110 J= UMAX 
GQ(JJC^,1) = GQ(J,K>N,1)/AA(J,K) 

GQ(J4C,N^) = GQ(J,KJSI,2)/AA(J,K) 

110 CONTINUE 


DO 100 N=U 

doiook=ucmax 

DO 100 J= UMAX 

GQN(J,K,N,1) = GQN(JJC^,1)/AAN(JJC) 
GQN(J,K t N^) = GQN(J,K,N > 2)/AAN(J ,K) 
100 CONTINUE 

RETURN 

END 


C*************** SCALAR TRIDIAGONAL ****** , •'* l •'*♦****♦****''" ,, *********** 


SUBROUTINE TRIV(JDIM,KDIMJLJU,KL,KUAB.C,F) 


DIMENSION A(JDIM,KDIM) 3(JDIM,KDIM),C(JDIM4Q)IM) 
DIMENSION X(JDIM3DIM)3(JDIM > KDIM) 


C 


DO 10 J=JL JU 

X(J,KL)=C(J,KL)/B(J t KL) 

F(J,KL)=F(J 1 KL)/B(JXL) 
10 CONTINUE 


KLP1 = KL +1 


DO 1 I=KLP1,KU 
DO 20 J=JL,JU 
Z=17(B(JJ)-A(J,I)*X(JJ-1)) 
X(JJ)=C(J,I)*Z 

F(JJMF(J,D-A(JJ)*F(J,I-1))*Z 
20 CONTINUE 

1 CONTINUE 
C 


KUPKL=KU+KL 


DO 2 11=KLP1JCU 
I=KUPKL-I1 


DO 30 J=JLJU 
F(JJ)=F(JJ>X(JJ)*F(JJ+1) 
30 CONTINUE 

2 CONTINUE 
C 

RETURN 

END 
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£***********+********************************************************** 
c *************** SCALAR PERIODIC TRIDIAGONAL *********************** 

0 *********^*********************************************************** 
SUBROUTINE VTRIP(JDIM,KDIM J 1 J2.K1 ,K2, A3,C 

c 

DIMENSION A(JDIMJCDIM)3(JDIM,KDIM),C(JDIMJCDIM)J(J D I M ,KD IM ), 
> QD(JDIM^DIM),S(JDIM > KDIM)J=ND(KDIM) 

JA = J1 + 1 

c 

C FORWARD ELIMINATION SWEEP 

C 

DO 1 K = K1.K2 
FND(K) = F(J2,K) 

QD(J1,K) = -C(J1,K)/B(J1,K) 

F(J1,K) = F(J1,K)/B(J1,K) 

S(J1,K) = - A(J1,K)/B(J1,K) 

1 CONTINUE 
C 

DO 10 J=JAJ2 
DO 2 K = K1.K2 

P =17( B(JJC) + A(J,K)*QD(J- 1 ,K)) 

QD(J,K) = - C(J,K)*P 

F(J,K) = ( F(JJC) - A(J,K)*F(J-UC))*P 

S(J,K) = - A(J,K)*S(J-1,K)*P 

2 CONTINUE 

10 CONTINUE 
C 

C BACKWARD PASS 

JJ = J1 + J2 
DO 3 K = K1,K2 
QD(J2JC) = 0. 

S(J2,K) = 1. 

3 CONTINUE 
C 

DO 11 I=JAJ2 
J = JJ-I 

nn A V: — V 1 V') 

S(JJC) = S(J,K) + QD(JJC)*S(J+1,K) 

QD(J,K) = F(J,K) + QD<J ( K)*QD(J+1,K) 

4 CONTINUE 

11 CONTINUE 
DO 5 K = K1.K2 

F(J2,K) = ( FND(K) - C(J2,K)*QD(J1,K) - A(J2,K)*QD(J2-1 JC))/ 

1 ( C(J2,K)*S(J1,K) + A(J2,K)*S(J2-1JC) +B(J2JC)) 

5 CONTINUE 
C 

C BACKWARD ELIMINATION PASS 

DO 12 I=JAJ2 
J = JJ -I 

DO 6 K = K 1 ,K2 

F(JJC) = F(J2,K)*S(JJC) + QD(J»K) 

6 CONTINUE 

12 CONTINUE 
RETURN 
END 
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c************************************************************** 


C********>m>******metoICS subroutine************************* 


C************************************************************** 


SUBROUTINE XYMETS(JMAX,KMAX,X,Y,XY,XYJ) 

C 

PARAMETER(JDIM= 1 20,KDIM=80,MAXJ= 1 2 1 ,MAXK=8 1 ) 
DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM),XY(JDIM,KDIM,4) 
DIMENSION XYJ(JDIMJCDIM), WORKX(JDIM,KDIM), 

$ WORKY (JDIM,KDIM) 

LOGICAL PERIODIC, SHARP, CUSP 


READ IN X, Y COORDINATES 


** FOR MAT FOR ONE-ELEMENT GEOM** 

USE FORT. 10 FOR PLOT3D 

OPEN(12) 

READ(12,*) JDUM.KDUM, AL 
READ(12,100) ((X(J,K), J=1 JDUM), K=1,KDUM) 

READ(12,100) ((Y(JJC), J=1JDUM), K=1JCDUM) 

REWIND 12 
100 FORMAT(10E19.12) 

C ** FORMAT FOR TWO-ELEMENTS GEOM** 

C OPEN(12) 

C READ(12 ,*) JDUM.KDUM 

C READ(12,*)((X(JJC),J=1JDUM),K=1JCDUM), 

C % ((Y(JJC), J=1JDUM), K=1,KDUM) 

C REWIND 12 

C 

C CHANGE DIRECTION OFT TO CLOCKWISE FOR BB MODEL 

C 

DO K = 1, KMAX 
DO J = 2, JMAX 
WORKX(J,K) = X(J,K)/AL 
WORKY(JJC) = Y(JJC)/AL 
ENDDO 

X(UQ=X(1JC)/AL 
Y(1,K)=Y(1 JC)/AL 
ENDDO 
C 

DO K = 1, KMAX 
DO J = 2, JMAX 
X(J,K) = WORKX(JMAX+2-J JC) 

Y(J,K) = WORKY (JMAX+2-J JC) 

ENDDO 

ENDDO 

C 

JLOW=l 
JUP=JMAX 
C JLOW=2 

C JUP=JMAX-1 

C XY4 =XXI,XY3 = YXI 

C 
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DO 11 K=1,KMAX 
DO 10 J=JLOW j JUP 
JP1 = J + 1 
JM1 = J - 1 

EF(J .EQ. 1) JM1 = JUP 
IF(J .EQ. JUP) JP1 = JLOW 
XY(J,K,4) = ( X(JP1,K) - X(JM1,K))*.5 
XY(J,K,3) = ( Y(JP1,K) - Y(JM1JC))*.5 

10 CONTINUE 
C 

JM=JMAX-1 

C IF(.NOT.PERIODIC)THEN 

C XY(1,K,4) = ,5*( -3.*X(1,K) +4.*X(2,K) - X(3,K)) 

C XY(1JC,3) - .5*( -3.*Y(UC) +4.*Y(2,K) - Y(3,K)) 

C XY(JMAXJC,4) = ( 3*X(JMAX,K) -4*X(JM,K) + X(JM-1JC))*.5 

C XY(JMAX,K,3) = ( 3* Y(JMAX JC) -4.*Y(JM,K) + Y(JM-1JC))*.5 

C ENDIF 

C 

11 CONTINUE 
C 

C XY2 = XETA, XY1 = YETA 

C 

DO 21 J= UMAX 
DO 20 K=2,KMAX-1 
XY(J,K,2) = ( X(J,K+1) - X(JJC-1))*.5 
XY(J,K,1) = ( Y(J,K+1) - Y(J,K-1))*.5 

20 CONTINUE 

XY(J,1,2) = ( -3.*X(J,1) 44.*X(J,2) - X(J,3))*.5 
XY(J,1,1) = ( -3*Y(J.l) +4.*Y(J,2) -Y(J,3))*.5 
KM=KMAX-1 

XY(JJCMAX^) = ( 3.*X(JJCMAX) -4*X(JJCM) + X(J,KM-1))*.5 
XY(JJCMAX,1) = ( 3.*Y(J,KMAX) -4.*Y(J,KM) + Y(J,KM-1))*.5 

21 CONTINUE 
C 

C FOR PERIODIC GRIDS WITH SHARP OR CUSP TRAILING EDGES USE 
C FIRST ORDER DERIVATIVE FOR ETA TERMS , SECOND ORDER 
C SOMETIMES LEADS TO NEGATIVE JACOBIANS 

PERIODIC = .TRUE. 

SHARP = .TRUE. 

CUSP = .FALSE. 

IF( PERIODIC )THEN 
J= 1 

IF( SHARP .OR. CUSP)THEN 
XY(J.U)= -X(J,1) + X(J,2) 

XY(J,1,1) = -Y(J,1) + Y(J^) 

ENDIF 

ENDIF 

C 

DO 30 K=1JCMAX 
DO 30 J= UMAX 

DINV = 1./ ( XY(J,K,4) * XY(J,K,1) - XY(J,K,3) * XY(J,K^) ) 

IF(DINV .LE. 0.) THEN 
PRINT*, JJC 

PRINT*, XY(JJC,1),XY(JJC,2), XY(J,K,3),XY(JJC,4) 

ENDIF 

IF(DINV LE. 0.0)WRITE(6 ,*)'JACOBIAN( 'J.K,' )= '.DINV 
XYJ(J,K) = DINV 
30 CONTINUE 
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DO 32 K=1 JCMAX 
DO 32 J= UMAX 
DINV = XYJ(JJC) 

XIX = XY1, XIY=XY2,ETAX=XY3, ETAY = XY4 

XY(J,K,1) = XY(J,K,1)*DINV 
XY(J,K,2) = - XY(J,K,2)*DINV 
XY(J,K,3) - - XY(J JC3)*DINV 
XY(J,K,4) = XY (J JC,4)*DINV 
32 CONTINUE 
C 

RETURN 

END 


c **,************ DE FTNmON OF MAJOR VARIABLES***************** 


AA WEIGHTS FOR VELOCITY GRADIENT COMPUTATION AT CELL CENTER 

AAN WEIGHTS FOR VELOCITY GRADIENT COMPUTATION AT GRID POINT 

AL CHORD LENGTH 

ALPI ANGLE OF ATTACK OF THE AIRFOIL 

BX LHS MATRIX ENTRIES IN THE XI DIRECTION 

BY LHS MATRIX ENTRIES IN THE ETA DIRECTION 

BWT INTERPOLATION RATIO OF SMIN LOCATION ON THE SURFACE 

CX LHS MATRIX ENTRIES IN THE XI DIRECTION 

CY LHS MATRIX ENTRIES IN THE ETA DIRECTION 

DAMP1 THE PRODUCT OF THE DAMPING FACTORS IN BB MODEL, D1 * D2 

DAMP1M THE PRODUCT OF THE DAMPING FACTORS IN BB MODEL, D1 * D2 

AT MID K POINT 

DAMP2 DAMPING FACTOR, F2(Y+) 

DS SPATIALLY VARIABLE TIME STEP, SCALED ON METRIC JACOBIAN 

DT TIME STEP 

DX LHS MATRIX ENTRIES IN THE XI DIRECTION 

DY LHS MATRIX ENTRIES IN THE ETA DIRECTION 

EDDY EDDY VISCOSITY IN ZETA 

FIRST LOGICAL: (T) THE SUBROUTINE IS CALLED FOR THE FIRST TIME 

FMU LAMINAR EDDY VISCOSITY 

FNU KINEMATIC VISCOSITY 

FSMACH FREE STREAM MACH NUMBER 

FX RHS MATRIX ENTRIES IN THE XI DIRECTION 

FY RHS MATRIX ENTRIES IN THE ETA DIRECTION 

GAMI GAMMA - 1 

GAMMA GAS CONSTANT (1.4) 

GQ VELOCITY GRADIENTS AT CELL CENTER 

GQN VELOCITY GRADIENTS AT GRID POINT 

B POINTER FOR SMIN 

ICST FLOW INDICATOR: 

0 - INITIAL RUN 1 - STEADY CASE 

2/3 - PITCHING UP 4 - OSCILLATING 

JDIM J DIMENSION OF THE GRID 

JMAX TOTAL NUMBER OF POINTS IN XI DIRECTION 

JTAIL1 FIRST SOLID BODY POINT 

JTAIL2 LAST SOLID BODY POINT 

KDIM K DIMENSION OF THE GRID 

KM AX TOTAL NUMBER OF POINTS IN ETA DIRECTION 

NNTT INDEX FOR ITERATIVE SOLUTION OF TURRE, NNIT = 1 FOR TIME 

ACCURATE SOLUTION 

PERIODIC PERIODIC(T) AND NON-PERIODIC(F) OPTIONS 
PI 4*ATAN(1.) 

PRESS PRESSURE (PLOT3D FORMAT) 

Q CONSERVATIVE VARIABLES (PLOT3D FORMAT) 

RE REYNOLDS NUMBER 

SMIN MINIMUM DISTANCE TO THE WALL 

S S PRODUCTION OF K OR TEMPORARY VARIABLE FOR BWT 

TINF INFINITY TEMPERATURE 

TURMU TURBULENT EDDY VISCOSITY 

TURRE TURBULENT REYNOLDS NUMBER 

U TANGENTIAL VELOCITY 

V NORMAL VELOCITY 

VORT VORTICITY 


X CARTESIAN X COORDINATES OF THE GRID 

XY METRIC TRANSFORMATIONS 

XY(JJC,1) = D XI /DX 
XY(JJK,2) = D XI / DY 
XY(J JC3) = D ETA / DX 
XY(J,K,4) = D ETA / DY 

XYJ JACOBIANS OF METRIC TRANSFORMATIONS 

Y CARTESIAN Y COORDINATES OF THE GRID 
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